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ABSTRACT 

We model the electromagnetic signatures of massive black hole binaries (MBHBs) with an associated gas 
component. The method comprises numerical simulations of relativistic binaries and gas coupled with calcu- 
lations of the physical properties of the emitting gas. We calculate the UV/X-ray and the Ha light curves and 
the Ha emission profiles. The simulations are carried out with a modified version of the parallel tree SPH 
code Gadget. The heating, cooling, and radiative processes are calculated for two different physical scenarios, 
where the gas is approximated as a black-body or a solar metallicity gas. The calculation for the solar metal- 
licity scenario is carried out with the photoionization code Cloudy. We focus on sub-parsec binaries which 
have not yet entered the gravitational radiation phase. The results from the first set of calculations, carried 
out for a coplanar binary and gas disk, suggest that there are pronounced outbursts in the X-ray light curve 
during pericentric passages. If such outbursts persist for a large fraction of the lifetime of the system, they can 
serve as an indicator of this type of binary. The predicted Ha emission line profiles may be used as a criterion 
for selection of MBHB candidates from existing archival data. The orbital period and mass ratio of a binary 
may be inferred after carefully monitoring the evolution of the Ha profiles of the candidates. The discovery 
of sub-parsec binaries is an important step in understanding of the merger rates of MBHBs and their evolution 
towards the detectable gravitational wave window. 

Subject headings: black hole physics — galaxies:nuclei — hydrodynamics — line:profiles — radiation mecha- 
nisms:general 



1. INTRODUCTION 

Massive black hole binaries (hereafter MBHBs) may form 
as a result of several p rocesses: through the b reakup of a 
supermassive protostar (Begelman & Rees1 fr978). as a meta- 
stable state in the evolution of a cluster of massive ob- 
jects dSaslaw et alJll974l:lBromm & Loebll2003f). and as a re- 
sult o f galactic mergers ([B egelman et al.lll980Hyaltaoia et alj 
[l989t IMilosavljevic & MerritJ 1200 It lYul 12001 ! The last 
mechanism is a major formation route for MBHBs and 
it relie s on hierarchical merger models of galaxy for- 
mation (Eaehnelt & Kauffmann 200l IVolonteriet all 120031) 
and significant dynamical evidence that the majority of 
galaxi es harbor massive black holes in their centers 
(e.g.. iKormendv & R ichstone 119951: iRichstone etail 119981: 
iPeterson & Wan del 2000). In this paper we investigate obser- 
vational signatures of MBHBs that result from galactic merg- 
ers. 

Theoretical results imply that the evolution of a MBHB 
proceeds in three stages. In the first stage, the binary sep- 
aration decreases through the process of dynamical friction 
when stars in the nuclear region are scattered by the binary. 
In the third stage, binary orbital angular momentum is 
efficiently dissipated via gravitational radiation emission and 
the binary proceeds to an inevitable merger. A possibility 
that in the second stage of evolution the binary merger stalls 
because of depletion of stars in the n uclear region, is still 
a top ic of a lively scient if ic debate dOuinlan & Hernquistl 
19971: ISigurdssorJ IT9981: IMilosavljevic & Merritt) 120011: 
Hems endorf etal.1 f2002k iMerrittl l2002t lYul l2002t karse'thl 



| 2003[ IMilosavljevic & Merritdl200l iPerets et alj 120061: Izieil 
2006). If this is indeed the case, binar ies spend most of their 
lifetime at separations of 0.01 — lpc (Begel man et al.lfT980l) 
and the majority of observed binaries should be found in 
this evolutionary stage. Hereafter, we refer to this, second 
evolutionary stage in the life of a binary, as to intermediate 
phase. 

A number of authors have suggested that interaction of a 
binary with a gaseous accretion disk may have si gnificant ef- 
fects on the inspiral and merger rates of MBHB (Ivanovet 



1999tlGould & Rixll2000l:lArmitage & Nataraianl|2002L 120051 
Escala et al.ll2004 120051: iKazantzidis et al J 120051 iDotti et al.1 
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2006at iMayer et alJ l2006aT) ~ It has been shown that mas 
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sive gas disks tend to form in centers of merging galaxie s 
dBarnes & Hernquistll 19921 119961: iMihos & Hernquistill996l) . 
In particular, a gas disk may expedite the merge r of a binary 
by dissipation of orbital angular momentum (lEscala et alJ 
12004 120051) . reducing the merger time to 10 6 - 10 7 yr. 

Observationally, MBHBs may be identified directly in cases 
when a black hole pair is spatially resolved and when the 
separation and kinematic parameters point to a close, bound 
pair. In practice, there is only one such object identified so 
far: observations of NGC 6240, an ultra luminous infrared 
galaxy, with the Chandra X-ray observatory revealed a merg- 
ing pair of X-ray activ e nuclei wi th a separation of 1 .4 kpc 
dKomossa et al.ll2003al) . Recently iHudson et all d2006l) also 
reported on the X-ray detection of a wide, proto supermassive 
binary black hole at the center of cluster Abell 400. Practical 
obstacles in this type of direct identification arise for several 
reasons. Firstly, fairly high spatial resolution and accuracy in 
position measurements are required to resolve a binary active 
nucleus. The spatial resolution required to resolve an interme- 
diate binary with orbital separation of lpc at the distance of 
^100 Mpc is about 2mas. With spatial resolutions currently 
achievable, it is easier to spot and resolve wide binaries, like 
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NGC 6240, whereas it has been suggested that MBHBs spend 
a majo r fraction of their life time as intermediate, hard bi- 
naries (Begelm an et al.l [l980). Because of their high spatial 
resolution VLBA and VLBI radio observatories have a large 
potential for discovery of close binaries with separations of 
0. 1 — 10 pc, especially when combi ned with ground-bas ed op- 
tical spectroscopic observations (Rodriguez et al. 2006|). 

The second most appealing piece of observational evidence 
for a MBHB would be a periodic manifestation of Keplerian 
motion. There are a handful of galaxies which exhibit pe- 
riodicities in their ligh t curves and thus can be considered 
as MBHB candidates dFan et al.lfl998t iRieger & Mannheim! 
2000; De Paolis et all 120021 12003b LSudou et al.l 12003b IXid 
2003b . A prominent example is the blazar OJ 287 which ex- 



120031) . These are expected to be an important class of 
sources for the Laser Interferom e ter Space Anten n a, LISA 
dBenderi [19981: iDanzmannl [l996t lHaehneil [l994t iHughesI 
2002b IWvithe & Loebl 12003b ISesanaTtal] |2004 120051: 
Rhook& Wvithe 2005). Since it is not clear which model 



hibits outburst activity in the optical light curve with a pe- 
riod close to 12 years. The latest outburst of OJ 287, ex- 
pected at the end of 2006, has not yet happened, putting to 
the test a number of models along with the binary hypothesis 
dSillanpaa et al.lll988blLehto & Valtonen1ll996blVartaoja et alJ 
120001) . 

Another phenomenon associated with the existence of 
a MBHB is the change in the orientation of the black 
hole spin axis and its precession about the total angu- 
lar moment um axis, caused by its interaction with another 
black hole dMerritt & Ekersl 120021: IZier & Biermamil 1200 ll 
2002). This effect is one possible explanation for helical 
radio-jets, X-shaped , and S-shaped radio-jets, observed in 
some galaxies (e.g.. | Hunstead et ail 1984b iLeahv & Wi lliams 
1984b IParma etaO 1985b IRoos et al. 119931 IWang et alJl2003b 
Liu! 12004}). though alternative explanations are available 
dCapetti et alj|2002t) . In cases when the orbit of the less mas- 
sive, inspiralling black hole becomes coplanar with an accre- 
tion disk of the primary, fueling of jets may be interrupted 
and re-started dLiull2004T) . This could give rise to a second 
pair of jets aligned along the same axis with the initial pair, 
as seen in double-double r adio galaxies dSchoenmakers et al] 
l2000blLiu. Wu. & Caoll2003l) . 

Massive binary black holes may also influence the kine- 
matics of the broad line region of the active galactic nu- 
cleus (AGN) which coul d be detected in per i odic variation s 
of broad emiss ion lines ([Begelman et al.ll980llGaskei]|1983l) . 
Gaskell (1996) suggested that double peaked emission lines 
originate in binary broad line regions. In this scenario peri- 
odic velocity shifts of each peak in the line profile should re- 
flect the orbital motion of a single broad line region associated 
with one of the black holes. After extensive spectroscopic 
monitoring of MBHB candidates spanning two decades the 
predicted long term periodic variability was not confirmed, 
leaving little room for a bin a ry broad line region scenari o 
dHalpern & Filippenkd I1988L 1 1992b lEracleous et ail 119971) . 
lEracleous et alJ (119971) noted that perturbation of a single 
broad-line region by a second MBH is still an open possibility, 
and in fact the only good interpretation for some objects. An- 
other version of this hypothesis in which the secondary black 
ho le perturbs the broa d line region of the primary was offered 
bv lTorres et alJ (l2003l) as an explanation of the variab le Fe Ka 
line in 3C 273. Etherin gton & Ma cieiewskil d2006l) suggest 
that characteristic spiral arm morphology can arise within the 
inner 1 kpc of a gas disk as a consequence of perturbations 
in a binary potential, and that in the absence of more direct 
techniques, this signature may be used for a MBHB search. 

Another, soon-to-be-opened observational window 
should allow detection of gravita t ional waves from MB- 
HBs dThorne & Braginskiil Tl976fc ICentrellal l200l iBaker! 



should be used to describe the optical data (binary orbital 
motion or disk and jet precession), the masses of current 
candidate binary systems are highly uncertain. The future 
data from gravitational radiation observations, combined with 
information on optical counterparts, should play an important 
role in determining the most like ly scenario dKomossal 
2003b; Armitage & Natarajan 120051) . In a broader sense, 
once enough information about the statistical properties of 
such binaries is available, it should be possible to constrain 
the merger rates and merger h istories of massive and super- 
massive black hole binaries dHaehneltl 1 1994b iMenou et alJ 
2001; Huehesl l2002h . For this reason a large scientific 
effort has been recently directed towards understanding 
the electromagnetic signatures in different pre-merger and 
post-merger phases of MBHBs (Milosavl jevic & Ph innev 
12005b iKocsis et all 120061: iDotti et al.l l2006bl) . The work we 
present in this paper fits within this broader context. 

In $2] we present our method for calculation of emission 
signatures of MBHBs. The results are described in $3] the 
discussion and future prospects in [jU and in |5]we give our 
conclusions. 

2. NUMERICAL SIMULATIONS 

We have carried out smoothed particle hydrodynamical 
(SPH) simulations of the binary and a gaseous component, 
and we have characterized the physical properties of the gas 
by calculating heating, cooling, and radiative processes as an 
integral part of simulations. Based on these results we have 
calculated the accretion-powered continuum and Ha light 
curves, as well as the Ha emission line profiles emerging 
from the inner parts of a gas disk on a scale of < 0. 1 pc. 

For the first time we attempt to model both gas and binary 
dynamics in the nuclear region on sub-parsec scales in order 
to bridge the gap between large scale hydrodynamic simu- 
lations of merging galaxies and numerical relativity simula- 
tions of the late, pre-coalescence evolutionary phase in the 
life of a binary. The binary-gas evolution in simulations of 
merging g alaxies has been followed down to separations of 
> 0. 1 pc (IDotti et aTll2006c|) . The simulations that reach even 
smaller separations are computationally expensive, mostly be- 
cause they need to be carried out at high resolution in force 
(or equivalently time). Additionally, if the nuclear processes 
followed include accretion, then a high resolution in mass is 
desired as well. These requirements are combined with the 
wide spatial and temporal dynamic range set by the nature of 
the problem itself. For the above reasons, in our simulations 
we focus on the length scale occupied by a nuclear accretion 
disk and do not consider the disk structure on l arger scales. 

We have used a mo dified version of Gadget dSpringel et alJ 
12001b ISpringell 120051) to carry out the MBHB simulations. 
Gadget is a code for collisionless and gas-dynamical cosmo- 
logical simulations. It evolves self-gravitating collisionless 
fluids with a tree A^-body approach, and collisional gas by 
SPH. Gadget was not originally intended to carry out rel- 
ativistic calculations. For this reason we have performed 
several modifications to the code in order to treat the two 
massive black holes relativistically. We have introduced 
the black hole particles as collisionless massive particles 
with pseudo-Newtonian potentials. The calculations with a 
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pseudo-Newtonian potential account for the decay of a black 
h ole b inary orbit through emission of gravitational radiation 
( £|2.21 >. The gravitational drag from gas particles on the two 
black holes is accounted for through a contribution from ev- 
ery gas particle to each black hole's acceleration. 

We set up the two Schwarzschild black holes as sink par- 
ticles and model the accretion rate of gas and resulting ac- 
cretion luminosity. The approach in which a black hole is 
treated as a sink particle ha s been introduced before in sim- 
ulations of star formation dBate et al.l 120031: iLi et alj 120051) 
and a ccretion onto supermassive black holes dSpringel et alj 
2005}). In our model, a particle is considered to contribute 
to the unresolved accretion flow near a black hole once it 
crosses the hole's accretion radius, Race- Race is set to a 
fiducial value of 20 r g , comparable to the smallest length 
scale resolved in the simulation (i.e., the smoothing length 
of particles in a region of high density, h m i ~ 10r g , where 
r g = GM BH /c 2 = 1.48 x 10 13 M 8 cm, and M BH = 10 8 M 8 M 
is the mass of a black hole). Once the gas particles find them- 
selves within Race their dynamics is not temporally resolved. 
This choice allows to avoid a simulation hang-up when the in- 
tegration time step of such particles becomes very small. We 
introduce a numerical approximation in the calculation of the 
accretion rate in order to account for the fact that the unre- 
solved accretion flow of particles inside R ace radius is noisy 
and that their dynamics cannot be described by the model of 
a steady accretion disk. Because of the short dynamical time 
associated with the innermost region of the accretion disk the 
particles there respond quickly to the variable potential of the 
binary which causes the streams of gas to collide and intersect 
with each other. From the total amount of gas that crosses R acc 
in given time dt, only a fraction of the gas will be accreted, 
while the rest will gain enough momentum to avoid accretion 
onto the black hole. The accretion rate of a noisy flow of par- 
ticles across some radius r is in general case proportional to 
the surface area, M °< r 2 , assuming the planar and uniform 
distribution of particles within r. We calculate from the simu- 
lations the rate with which particles cross R acc and extrapolate 
from it the accretion rate of particles that would have crossed 
the Schwarzschild radius as M(2r g ) w M(R acc ) (R acc /2r g )~ 2 . 
The estimate obtained in such way is conservative as we do 
not take into account the gravitational focusing by the black 
holes, which would increase the flow of particles towards the 
black holes and thus, their accretion rates. In §|4] we discuss 
the effect of our assumptions about the accretion model on the 
observational signatures. 

2.1. Viscosity 

The viscosity is parametrized in numerical calculations in 
order to account for the transfer of momentum and angu- 
lar momentum without detailed assumptions about the na- 
ture of the phenomenon. One widely used parametrization 
is based on the Shakura-Sunyaev model ( Shakura ~& Sunvaevl 
1 1973b . where the a-par ameter takes values less than 1.0, 
and typically a < 0.1. Starling et al.l (120041) for example 
find 0.01 < a < 0.03, for accretion disks in AGN, based 
on the observed continuum variability. Finding a suitable 
parametric description for the viscosity in SPH calculations 
which would faithfully describe the angular momentum trans- 
fer in different physical regimes has been a major issue 
over many ye ars dHernquist &K atz 1989; Monaghan 1991 
Balsaral [1991 ISteinmetzl [l996t iMonaghanl Il997t ISpringell 
20051) . A notable improvement in the treatment of viscos- 



serving the en tropy of a simulated system (e.g., in Gadget-2; 
Springel 2005). Although we have us ed the version of Gadget 
released earlier (Springel et al. 2001), for the purposes of our 
calculations we adopted the viscosity prescription given in the 
new version of the code, Gadget-2. The physical properties in 
our simulations that are directly affected by this choice are the 
accretion rate of gas onto black holes and the internal energy 
of the gas — both of which play very important roles in the 
observational appearance of the system. 

We derive the value of the artificial viscosity parameter used 
in the code, OCcadget m such wa Y mat it corresponds to the 
Shakura-Sunyaev viscosity parameter a < 1. We note that 
this relation is only a convenient way of describing the viscos- 
ity numerically and should not be interpreted as a statement 
about the physical viscosity in the disk. In SPH simulations, 
the viscosity acts as an excess pressure which is assigned 
to particles in the equation of motion, P « Ucadget (c s CO — 
jC0 2 )p, where c s is the speed of sound and CO is the rel- 
ative velocity of the two approaching particles. We com- 
pare this to the Shakura-Sunyaev tangential stress component 
P w —ac 2 p, and obtain a local relation between the parame- 
ters describing a single particle interaction 



CO 



3oj2 

2 c? 



(1) 



The resulting artificial viscosity as described by equation (fl} 
does not explicitly depend on the density of gas particles. 
However, in simulations a major increase in internal energy 
due to viscous dissipation occurs in high density regions 
which are conducive to high rate of particle interactions (i.e., 
collisions). The implication is that the relation connecting the 
two viscosity parameters depends on the density of gas parti- 
cles in a simulation, namely, 



a « OCcadget 



CO 

Cs 



3 0)2 

2 c 2 



(2) 



where p = p(r) is the density distribution of gas particles and 
Po = p{ r in) is the value of the density at the inner radius of 
the accretion disk. Assuming that the largest value that CO 
can take is of the order of particle orbital velocities one can 
find a value of Otcadget which corresponds to some value of 
the Shakura-Sunyaev a. From equation (0) we find that in 
order to have 0.01 < a < 0.1, at radial distances £, ~ 100 
from a black hole, in a disk where c s ^100 kms~', we have 
to choose CtQad^et ~ 10~ 6 . Hereafter, we use ^ as the radius 
in units of gravitational radii (£, — r/r g ). 

2.2. Elliptical Orbits in the Paczynsky-Wiita Potential and 
Gravitational Wave Emission 

We choose a pseudo-Ne wtonian Paczynsky-Wiita potential 
dPaczynsky & W iita 1980) to describe the potentials of mas- 
sive black holes in our calculations 



-GMbh 
(r~r s ) 



(3) 



ity was offered by parametrization that perform better in con- 



where r s — 2GMbh / c 2 ■ The two black holes can initially be 
placed on circular or elliptical orbits. The properties of cir- 
cular orbits in the Paczyn sky-Wiita potential are described in 
iPaczynsky & Wiital (ll980). In addition, we derive expressions 
for the orbital velocity of a test particle on an elliptical orbit 
in the Paczynsky-Wiita potential in order to assign initial ve- 
locities consistent with the potential. We adopt the reduced 
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two-body problem parametrization where M = m\ +rri2 and 
jj. =mim2/{mi+m2), where the energy and angular momen- 
tum of a system are calculated with respect to its center of 
mass (;«i and ni2 are masses of the two black holes). We start 
with the requirement that in the Paczynsky-Wiita potential, 
the total energy and orbital angular momentum at the peri- 
center and apocenter of the orbit are equal. This leads to the 
following conditions: 

1 1 (4, 



vl - vl = 2GM 



r a - r s 



r p V p = r a V a , (5) 

where v p , V a , r p , and r a are the velocity and radial distance of 
the reduced mass body, ju, at the pericenter and apocenter, re- 
spectively. From equations <j4j and © we derive the velocity 
at the pericenter as 

/ 2GM Vl 



"i- \ , 6 \ • ( 6 ) 

y r a + r p (r p - r s ) [r a - r s ) 

From the expression for the orbital angular momentum, 
L/n = r p v p = r 2 9 (where r is the radial distance of the re- 
duced mass from the center of mass of the system, and 9 is its 
true anomaly), we find the angular azimuthal velocity. 



9 



1 2GM 



r r 
' a' p 



r a + r p {r p -r s )(r a -r s ) 



(7) 



We derive the radial velocity component from the equation of 
an ellipse, r = a(l — e 2 )/(l — ecos 9), 



dr 
di 



dr 



-9 



e sin 9 r 9 



(8) 



d9~ a{\-e 2 ) 

Using the equation (|6]l and the condition E(r p ,v p ) = E(r,v) 
we find the vis viva expression for an elliptical orbit in the 
Paczynsky-Wiita potential to be 

V = (2GM)'/ 2 x 

11 1 1 1 \ 1/2 

(9) 



<p-r s (r„/r p ) 2 ~l r a ~r s {r p /r a ) 2 -\ 
As the binary orbits in Paczynsky-Wiita potential it loses a 
fraction of its energy and orbital angular momentum through 
gravitational wave emission. We assume that as the binary 
contracts by emission of gravitational waves, the two black 
holes remain on elliptical orbits. This assumption is justified 
in our simulations where the total mass of the gas in a simula- 
tion is much smaller than the mass of a binary. Consequently, 
the perturbations of the elliptical orbits of black holes by the 
underlying g as distribution are small. 

Following lLandau & Lifshitz] (1 19751) we calculate the rate 
of the energy loss to gravitational radiation using the 
quadrupole approximation 

dE _32C^W 

dt 5c 5 ' 1 j 

where £2(r) is the orbital frequency in the Paczynsky-Wiita 
potential. From this equation it is possible to obtain the radi- 
ation r eaction acceleration at the center of mass frame of the 
binary (Lee & Kluzniak 1999, and references therein) 



1 



1 dE v 



(11) 



fi dt v 2 

chosen so as to preserve the symmetry of forces acting on each 
member of the system. This requirement ensures that the total 
momentum of the system is conserved while the individual 
components recoil due to the emission of gravitational waves. 



2.3. Heating and Cooling of the Gas 

We have investigated two different scenarios for radia- 
tive heating and cooling, assuming the emitting gas can 
be described as a black-body (BB) or a solar metallic- 
ity gas. The option for radiative heating and cooling of 
hydrog en-helium gas is also implemented in the code, fol - 
lowing jThoul & Weinberg] dl996l) . iBlackl dl981l) . ICeiil dl992l) . 
and iKatz et alj (1996), but was not used in the runs reported 
here. The cooling scenarios are mutually exclusive and by 
comparing them one can study the importance and impact of 
various cooling mechanisms in our simulations. In our sim- 
ulations, the gas is heated by shocks and "external" sources 
of illumination. We assume that the sources of illumination 
are powered by accretion onto the massive black holes. After 
accretion on either of the two black holes becomes signifi- 
cant, UV and soft X-ray radiation emitted from the innermost 
portion of the accretion flow photo-ionizes the gas. The bolo- 
metric luminosity of an accretion-powered source of ionizing 
radiation can be written as: 



GM B h m 1 . 2 

Lace = 1 ^ = r]Nm p c z 

2r s 2 



(12) 



where 77 = 0.01 is the assumed radiative efficiency, m and 
are the mass and particle accretion rates, respectively, and m p 
is the mass of a gas particle. 

1 . Black-body case. 

In this scenario, the illuminated gas absorbs and emits 
the radiation as a black-body. Each gas particle is there- 
fore treated as an optically thick cell of gas, in local 
thermal equilibrium within a radius equal to its smooth- 
ing length, h sm [. A fraction of the incident energy ab- 
sorbed by the gas depends on a coverage factor of the 
gas £, i.e., the fraction of solid angle covered by gas, 
as seen from the position of a source, and the distance 
of an absorbing gas cell from the source, r. From our 
simulations we find £ ~ 10~ 4 , corresponding to a frac- 
tion of particles exposed to illumination from an ex- 
ternal source at any given time. We do not take into 
account radiation pressure, i.e., the momentum that is 
transferred to the gas particles by the incident radiation. 
We write the expression for the heating rate of the gas 
per unit volume as 



bb 



2 m„ 



erg cm 3 s 1 . (13) 



Here we assume that Lx is the portion of the accretion 
luminosity carried by ionizing photons (Lx /L acc w 0. 1). 
The black-body cooling rate per unit volume is evalu- 
ated from the Stefan-Boltzmann equation and it can be 
written as 



A bb = 4nh 2 ml oT A — erg cm 3 s 1 , 
trip 

where a is the Stefan-Boltzmann constant. 



(14) 



2. Solar metallicity gas. 

Spectroscopic observations of AGN show a wide range 
of metal absorption and emission lines, implying that 
most AGN systems are metal-rich. Additionally, met- 
als are important coolants under the physical conditions 
typically found in AGN broad line regions. Therefore, 
in order to correctly model and describe these systems 
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it is necessary to also include metals as constituents 
of the gas. However, detailed calculations of radia- 
tive transfer, coupled with parallel hydrodynamic sim- 
ulations present a large computational challenge and 
are currently not possible. Including metals in heat- 
ing and cooling calculations results in a large network 
of equations, especially if non-LTE processes are in- 
cluded. The solution of the heating and cooling equa- 
tions dominates the CPU usage in a parallel hydrody- 
namics code to the extent that the calculation becomes 
prohibitively expensive. 

Here we describe the approximate method we used for 
calculation of heating and cooling of the gas with met- 
als. We have constructed a number of cooling maps 
using the photoionization code Cloudy dFerland et alJ 
1998). Since the maps are used in computationally 
expensive numerical calculations, it is impractical to 
call Cloudy every time a cooling rate is needed. In- 
stead, we pre-calculate a grid of cooling maps over a 
wide range of parameters, where the range of param- 
eters has been determined from preliminary calcula- 
tions without cooling. The parameter grid is read in 
by the simulation code, and radiative heating and cool- 
ing rates are linearly interpolated from the existing grid 
points. The cooling maps are calculated in the param- 
eter space of density and temperature of the gas and 
intensity of ionizing radiation. The range of parame- 
ter values for which the maps were computed is as fol- 
lows: 10 9 cm- 3 < n < 10 19 cm- 3 , 2000K < T < 10 8 K, 
and erg cm~ 2 s _1 <J < 10 17 erg cm _2 s _1 . Because 
one of the assumptions in Cloudy is that the electrons 
are non-relativistic, the present range of its validity ex- 
tends to temperatures below roug hly 10 9 K. Gas at 
higher temperatures than this is commonly encountered 
in our simulations; in such cases we calculate heating 
and cooling rates by linearly extrapolating the grid val- 
ues. In the SPH simulations we set the lower threshold 
for the gas temperature to 100 K and assign no upper 
threshold. The calculation in Cloudy is explicitly one 
dimensional, with results depending only on the depth 
coordinate. Propagation of line photons is conducted 
using the escape probability approximation. In this ap- 
proach the escape probability for a photon is calculated 
under the assumption that the difference between the 
mean intensity averaged over the line and the source 
function at some location in the e mitter is due to pho- 
tons leaking away from the region (Ferland 2003). Al- 
though this approximation is used in many state-of-the- 
art plasma codes because of its computational facility, 
it is not known to what extent it affects predictions of 
the line intensities. 

The density and temperature of the gas are evaluated in 
the SPH simulation after every time step. The spectral 
energy distribution (SED) of the ionizing radiation is 
taken to be a power-law with an index a = 1, extend- 
ing from Vi = 1.36 eV to V2 = 100 keV, a distribution 
commonly assumed for the central engines of AGN. 

J v =J (jpj ergcm-V'Hz- 1 , (15) 

where Vo is some arbitrarily chosen frequency and Jq = 
J(Vq) is the normalization of a SED that can be directly 
evaluated from the luminosity of an accretion source, 



eq. (fT2l . as follows 

Face = j^r = Jv dv » 1 1 .2 7 Vo , (16) 

where F acc is the total flux incident on the gas cell lo- 
cated at the distance r from the ionizing source. In ad- 
dition to the cooling processes included in Cloudy we 
also consider Compton cooling from a thermal distribu- 
tion of non-relativistic electrons, which may become a 
significant coolant for the hottest and densest gas in the 
nuclear region. 

Ac omp = 4-48x10-^) ( Ip £ Ff ) 

/ F a cc \ -3-1 /n\ 

x TTTTo T ergem-s 1 , (17) 

\10 lu ergcm z s l J 

where n e is the electron number density. When infor- 
mation about the cooling rate per unit volume, A* = 
^Cloudy k + h-Comv k > is obtained for every gas particle k 
in the simulation, their internal energy is modified ac- 
cordingly. Thus, the heating and cooling are directly 
coupled with the evolution of the physical conditions 
in the gas. The total cooling rate from the gas (i.e., the 
bolometric luminosity) is obtained by summing over all 
particles, "" A# (m p /pk)- We do not assume that the 
gas is in thermal equilibrium, and consequently cooling 
is not necessarily balanced by heating. The integration 
scheme used in implementation of cooling for both the 
black-body and solar metallicity models is explicit. The 
cooling time for individual gas particles in our simu- 
lations is typically longer than the time step assigned 
based on their acceleration, dt <* 1/ a p , and the Courant 
condition. Even when the cooling time is very short, 
dt < 10 2 s < 

tcooh guarantying the numerical stability 
of the integration. In all of the calculations the default 
abundance pattern is assumed to be solar. The atomic 
species considered by Cloudy range from hydrogen to 
zinc. 

Since Cloudy includes a detailed treatment of radiation 
processes, we have also used it to calculate the Ha line 
intensity and optical depth, the electron scattering opti- 
cal depth, and the neutral hydrogen column density for 
each gas cell. As we will describe in the next section 
this information allowed us to relax some of the impor- 
tant assumptions commonly used to calculate the Ha 
light curves and emission line profiles. 

2.4. Calculation of Ha light curves and emission line profiles 

We chose the Balmer series Ha line (X rest = 6563 A) to de- 
scribe the emission signatures of a MBHB. The broad Ha 
line is thought to reflect the kinematics of t he gas in the 
broad line region of AGN (ISulentic et a l. 2000), and it is typ- 
ically found to emerge from a wider range of radii in the 
broad line regio n than the UV and X-ray band broad emis- 
sion lines (e.g.. lEracleous et al.l Il996t IWandel et al] 1 19991 : 
ISulentic et af] l2000). Also, the Ha line is the most promi- 
nent broad line in the optical spectrum and the least contam- 
inated by the neighboring narrow emission lines. It is mostly 
powered by illumination of a di sk by an emission source 
(ICollin-So uffrin & Dumont 1989). The emission sources in 
our simulations are the innermost regions of the accretion flow 
around the two massive black holes that give rise to ionizing 
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UV and soft X-ray radiation (Shakura & Sunvaevll973l) . The 
cool gas is optically thick to the incident UV/X-ray photons 
which can escape the disk only after they have been converted 
to photons of lower energy. One portion of the absorbed 
UV/X-ray radiation is reproces sed by the gaseous disk and re- 
emitted as the broad Ha line (Shakura & Su nyaeviri973l see 
their discussion of broad recombination and resonance emis- 
sion lines). 

The most important consideration in calculating the emit- 
ted Ha light is the efficiency with which the gas repro- 
cesses the incident ionizing radiation and re-emits it in the 
Ha band. This efficiency is commonly characterized by the 
surface emissivity of the gas. The emissivity of the photoion- 
ized gas depends on numerous physical parameters. Locally, 
it depends on the physical properties of gas. Globally, its spa- 
tial distribution depends on the structure and morphology of 
the accretion disk. The Ha emissivity can be most accurately 
assessed with photoionization calculations, thus we used this 
approach in the case of a solar metallicity gas. 

In cases when photoionization calculations are not avail- 
able, the emissivity is often parametrized as a function of 
radius, e = £o^~^, where ^«3. This is justified by the 



1 s I 1 1 1 1 1 ' ' r 




10" 6 - 



j,Qr*l , i , i , . , i , . , 

2000 4000 6000 

Radius (r g ) 

FIG. 1 . — Ha emissivity of the gas as a function of radius, at 9.4 years for a 
model with solar metallicity gas. The emissivity of each gas cell plotted in the 
figure is weighted by the density of the gas at that position, so that comparison 
with the parametric emissivity model (plotted as a solid line with arbitrary 
normalization) can be made. The temperatures of gas particles are marked 
with color: red T < 10 4 K; yellow 10 4 K < T < 10 s K; green 10 6 K < T < 10 8 
K; blue 10 8 K < T < 10 10 K; and violet T > 10 10 K. This figure corresponds 
to the morphology of the disk plotted in panel 2 of Figure |7J 



photoi onization calculations of ICollin-Souffrin & Dumontl 
(1989), which suggest that the emitted Ha flux is approxi- 
mately proportional to the incident continuum flux for a wide 
range of values of the density and column density. We have 
compared the emissivity determined from our own photoion- 
ization calculations with the parametric model (Figure[T|) and 
found a qualitative agreement between the two. We also no- 
tice a significant amount of scatter in the emissivity values 
from the photoionization calculation, which is expected due 
to variations in surface density of the perturbed gas disk in 
our model. 

Once the surface emissivity of each gas cell in the sim- 
ulation is known, the Ha luminosity can be expressed as 

Lhoc x HtS* e k> where A: is a particle index. The implicit 
assumption, following from the above luminosity expression, 
is that the gas is optically thin to Ha photons. We find this as- 
sumption to be justified for a cold disk since its optical depth 
to Ha light is less than 10~ 4 . We show the optical depth of 
gas cells to the Ha photons as a function of the disk radius 
in Figure |2] When the gas is exposed to an intense ioniz- 
ing continuum, the optical depth to the Ha photons increases 
significantly because many hydrogen atoms are in the first ex- 
cited energy state. However, in our profile calculations we do 
not correct for the absorption of the Ha photons by the gas. 
This is because we are not able to correct consistently for ab- 
sorption of photons along the line of sight with the numerical 
code we have used. The effects of our assumptions on the Ha 
luminosity and emission line profiles are discussed in § !3.5l of 
the text. 

The Ha emission line profiles have been calculated 
taking into account the relativistic Doppler shift, v„/ M = 

Vrest \J 1 — fi 2 / (1 + P cos 0), and the gravitational redshift in 
the potential well of a Schwarzschild black hole, v b s = 

V rest yl —2/<i;, where V b s and v rest are the observed and rest 
frequencies of an emitted photon, j3 is the velocity of an emit- 
ter in units of speed of light, and 9 is an angle between the 
direction of motion of an emitter and the observer's line of 
sight. In the next section we summarize our choice of param- 
eters used in SPH simulations as well as the calculation of 
light curves and line profiles. 

2.5. Choice of Parameters 

In the calculations presented here, both black holes 
are treated as massive particles in the pseudo-Newtonian 
Paczynsky-Wiita potential. The masses of the primary and 
secondary black holes, together with the other parameters are 
shown in Table [T] In all simulations the binary mass ratio is 
chosen to be 0.1. 

Based on current theoretical and observational understand- 
ing of MBHBs it is not obvious what kind of orbits the bina- 
ries have in the intermediate phase of their evolution. The re- 
sulting accretion rate and observational signature of a MBHB 
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FIG. 2. — Left: Optical depth of gas cells to Ha photons as a function of radius at 9.4 years for a model with solar metallicity gas. The color legend is the same 
as in Figure □ This fi gure corresponds to the morphology of the disk as shown in panel 2 of Figure [JJ Right: Histogram showing the fraction of particles as a 
function of T(Ha). 

depend sensitively on the shape and orientation of its orbit 
with respect to the gas disk. The astrophysical significance 
of a predicted signature depends on the longevity of the par- 
ticular phase in the life of a binary and consequently on the 
probab ility that binaries will be detected i n that evolutionary 
phase (MacFadven & Milosavlievic 2006). The recent sim- 
ulations of iDotti et al.l (|2006a) imply that binaries which co- 
rotate with the gaseous disk circularize and loose the memory 
of their initial orbital eccentricity. On the other hand, counter- 
rotating binaries tend to preserve their original eccentricity 
down to w 1 pc, a limit set by the resolution of their simula- 
tion. In addition, simulations of gas disks with planets suggest 
that under some conditions tidally perturbed eccentric disks 
can exci te the eccentricity in a binary and t hat the two are 
coupled (Papaloizou. Nelson. & Masset 2001). Very few con- 
straints on the type of orbit can be placed from the observa- 
tions. An exception is the MBHB candidate OJ 287, for which 
it has been suggested that the binary orbit with and eccentric- 
ity of 0.7—0.8 and orbit al period close to 12 yr is coplanar 
with the accretion disk (ISillanpaa et al.l [19881: IValtaoia et al.l 
2000: lLiu&Wull2002l) . 

In the calculations presented here the primary and sec- 
ondary black holes are moving on elliptical orbits. At the 
beginning of a simulation the black holes start from the apoc- 
enters of their orbits (except in the model wh ere th e secondary 
starts with a true anomaly of 6 = 30°, see § I3.31 >. The initial 
values of the semi-major axis and eccentricity are 3007 r g and 
0.7, respectively. The period of the binary, calculated at the 
beginning of the simulation is 15.67 yr. 

The accretion disk is initially only associated with the pri- 
mary black hole. It is coplanar with the binary orbit and it 
extends from fy„ = 100 to % out = 2000 (% out translates to 0.01 
pc in physical units, a size expected for AGN nuclear accre- 
tion disks). The mass of the disk is 10 4 M Q and the number 
of gas particles in the disk is 10 5 . The mass of the thin disk 
is constrained by the requirement that its surface density be in 
the range that is commonly found in AGN broad line regions, 
where the broad Ha emission line is expected to originate. In 
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FIG. 3. — Evolution of the disk surface density and temperature profile 
before the secondary black hole is introduced into the simulation (IC2 run). 
The profiles correspond to progressively longer times as follows: years 
(solid, black line), 2.2 years (dotted, red line), 4.4 years (dashed, blue line), 
and 6.6 years (dash-dot, purple line) after the beginning of the simulation. 
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a Approximate times of pericentric passages. b Median value for the low and high density 
component. IC runs: median value for the unperturbed disk. c Estimated peak luminosity. 



addition, the disk mass is roug hly consistent with the results 
of the l arger scale simulat ions (Esc ala et~aT1l2005t iDotti et all 
l2006at iMayer et al.ll2006ah . where the mass of the gas disk 
enclosed within central 1 pc is of the order of the mass of the 
binary. Extrapolating from these results and assuming a con- 
stant surface density for simplicity, one obtains a disk mass of 
- 1O 4 M , for a disk size of 0.01 pc and a ~ 10 8 M© binary. 
In §|4] we discuss how the choice of the disk mass may affect 
the observational signatures. 

The initial rotational velocities of gas particles were calcu- 
lated based on the gravitational forces only while the effect 
of pressure was neglected in initially cold and unperturbed 
disk. We choose the initial temperature of the disk to be 
T = 2000 K everywhere. The initial surface density distribu- 
tion in the disk is £ °< t, , as set by the requirement that the 
disk has a uniform vortensity (defined as the ratio of vorticity 
to surface density). This criterion must be satisfied in order 
to construct the equilibrium model of an infinitesimally thin 
gaseous disk. The initial conditions have been chosen in such 
way that the disk remains near hydrostatic equilibrium before 
the secondary black hole plunges into it. Figure [3] shows the 
evolution of the disk surface density and temperature profile, 
from the initial conditions run with 100k particles (IC2) be- 
fore the secondary black hole is introduced. The gas in the 
IC2 run is assumed to have a solar metallicity. The tempera- 
ture of the disk settles around 1000 K in this period and there 
are no visible transients. No boundary conditions were ap- 
plied at the inner and outer edge of the disk in our simulations. 
Figure [3] shows a gradual evolution of the disk edges towards 
the regions of low density. This is a consequence of the fact 
that the disk has sharp edges, where particles are subjected to 
asymmetric pressure forces. We do not expect the diffusion 
of particles from inner edge towards the primary black hole 
to affect the results, since no particles are accreted during this 
initial period. At later times, in runs where the secondary 
interacts with the disk, the contribution to the accretion rate 
from a slow diffusion process should be negligible in compar- 
ison with the total accretion rate induced by the interactions. 

The viscosity parameter, radiative efficiency, and metallic- 
ity of the gas are kept constant, as described above; their val- 
ues can be found in Table Q] The softening parameter, <^ so f t , 
used in the expression for force in order to avoid numerical 
singularity, ranges between 0.1 and 0.2 (in units of Mg) for 
gas particles and is 4.0 for the massive black holes. 

We compute line profiles emerging from the gaseous disk 
under the assumption that the observer is located at a dis- 
tance d — > 00, in the positive xz-plane, at i = 30° to the z-axis. 



Changing the inclination, changes the values of the projected 
velocities (i.e., the overall width of the line profile) but has a 
relatively small effect on its shape otherwise. The azimuthal 
orientation of the observer is an additional parameter in the 
modeling of non-axisymmetric systems. The azimuthal orien- 
tation is measured in the orbital plane of a black hole binary, 
counter-clockwise with respect to the positive x-axis. In the 
calculations we present here, we adopt 0o = 0° as an arbitrary 
direction towards the observer. 



3. RESULTS 

We present and compare results of three different scenar- 
ios in which the evolution of a massive black hole binary and 
coplanar gaseous accretion disk is followed over two orbital 
periods of the binary (~ 30 years in total). The simulations 
span only two orbits because at the temporal resolution re- 
quired by the hydrodynamical calculation, we find it com- 
putationally expensive to follow a binary over many orbital 
periods. The calculation of the two orbits, with total of 10 s 
particles, required about 20 days of CPU time on 8 2.8GHz 
Intel Xeon processors per simulation. We find that 50-70% 
of the CPU time is consumed by the hydrodynamic calcula- 
tions, including the neighbor search, density determination, 
computation of hydro forces, and heating and cooling of the 
gas. A fraction of 20-30% of the CPU time is spent on com- 
putation of gravitational forces, including the tree walks and 
tree construction. We find 8 CPUs to be an optimal number 
of processors in terms of speed for calculations that employ 
10 s particles. Increasing the number of processors as 2" (i.e., 
to 16, 32, etc.) results in a loss of computational speed be- 
cause the CPU time consumption becomes dominated by the 
communication among the processors. 

In the first of three calculated scenarios the binary is co- 
rotating with the accretion disk which is assumed to absorb 
and emit radiation as a black-body. In the second and third 
scenario the binary is co-rotating and counter-rotating with 
respect to the disk, respectively, and the gas is assumed to 
have a solar metallicity. In the last two cases the heating and 
cooling of gas are followed with photoionization calculations. 
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FIG. 4. — Sequence of snapshots from the simulation showing the evolution of the binary and gas in the BB-model (projected into the plane of the binary orbit). 
The time stamp is marked on the top of each panel. The first snapshot shows the system shortly after the secondary plunges into the disk for the first time. The 
rotation of the binary and the disk is counter-clockwise. The color legend is the same as in Figure[T] Note that higher temperature particles are plotted over the 
lower temperature ones, with the result that some information is hidden. 



3.1. Co-Rotating Case With Black-Body Cooling (BB-model) 

The approximation in which emission from a geometri- 
cally thin accretion disk is represented by a modified black- 
body spectrum has been widely utilized since the work of 
Shakura & Sun vaevl dl973l) . We use this simple model as a 
control case for our simulations and we compare it with the 
photoionization models. 

Before the plunge of the secondary black hole into the disk 
the density in the quasi-steady disk ranges between 10 15 and 
10 12 cm at the innermost and outermost edges of the disk, 
and in the intermediate region it takes values according to 
the assumed power-law surface density. The temperature of 
the quiescent disk quickly adjusts to a two phase distribution. 
The higher density gas, which occupies the inner portion of 
the disk cools more efficiently (eq. [I4i i and reaches a temper- 
ature of ~ 10 2 K (Fig|4]l. The lower density gas occupies a 
distinct phase in temperature space at ~ 10 3 K. We are not 
able to follow the phase of steady accretion at low rates prior 
to the impact of the secondary. This is because our accretion 



disk is initially set with an inner edge at = 100, where vis- 
cous migration of gas in the radial direction is slow. For the 
set of parameters we have adopted, the theory of steady thin 
disks predicts an accretion rate of 10~ 4 — 10~ 3 M Q yr _ . The 
unperturbed disk is stable to gravitational collapse as deter- 
mined by the value of the Toomre parameter, Q = c s £1/ GE, 
where fi is the angular speed of the gas and E is its surface 
density. The innermost region of the gas disk has the lowest 
Toomre parameter, of order of a few but consistently larger 
than 1. The speed of sound, c s = {jp/p) 1 ^ 2 , also shows a 
characteristic two phase distribution with values lower than 
10 kms -1 across the disk, where 7=5/3 for a mono-atomic, 
ideal gas. The physical properties of the disk remain stable 
and do not change over a considerable period of simulation 
time, until the secondary black hole plunges into it. This is 
an indication that the disk has settled into a quasi-steady state 
and that any further physical changes can be attributed to the 
interactions with the binary. 
The tidal perturbation by the secondary causes the disk to 
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FIG. 5 . — Effective accretion rate onto the primary (solid, red line) and sec- 
ondary (dashed, blue line) black holes calculated from the BB-model. The 
accretion rate curves can be translated into UV/X-ray light curves, by as- 
suming that 1 M s yr~ 1 produces ~ 1 43 erg s~ 1 of UV/X-ray luminosity. The 
arrows mark the times of pericentric passages of the binary. 

become eccentric and develop a spiral pattern, even before the 
secondary plunges into it. The physical properties of the disk 
are significantly altered after the plunge, which occurs ~ 6 
years after the beginning of the simulation. The secondary 
forms a shock where the temperature persists at 10 12 K. While 
in this model the cooling of the gas is less efficient, there is no 
doubt that the shock has formed and that it trails behind the 
secondary as it travels through the disk. 

In our model, a source of illumination is assumed to turn on 
once the accretion onto a black hole begins. With both black 
holes accreting there are two sources illuminating and heat- 
ing the gas. As a result a portion of the gas is blown out of 
the plane of the disk and forms a hot "corona". The temper- 
ature of the gas is highest in the vicinity of the black holes, 
where it is in the range of 10 6 - 10 10 K. After the first pas- 
sage of the secondary the disk has formed a single trailing 
spiral which cools radiatively and adiabatically to T~ 10 3 K 
(FigHJi. The secondary black hole acquires a small, gravita- 
tionally bound corona of gas from the disk. After the second 
passage the secondary further heats and perturbs the disk gas 
which now reaches the scale of a parsec in the plane of the 
disk (recall the disk initial size 0.01 pc). The secondary black 
hole intersects the spiral arm of the disk and gravitationally 
captures a portion of it. For a while, the secondary moves on 
its orbit followed by the remainder of the spiral arm, which 
resembles a cometary tail. Since a large portion of the spi- 
ral arm, which retains a significant amount of colder gas in 
the disk, is disrupted after the second orbital passage, the ef- 
fective temperature of the gas becomes more dominated by 
the inner hot portion of the disk. The morphology of the sys- 
tem at the completion of the second binary orbit consists of 
a discernible denser disk component immersed into a spheri- 
cal halo of tenuous gas. The gas density now spans the range 



from 10 1U - 10 1 j cm J in the disk component and it extends 
down to 10 s cm~ 3 in the low density regions on the scale of 
~ 1 pc. The speed of sound is ~ 5 kms -1 and ~ 10 3 kms -1 
for the cold and hot gas phases, respectively (Table|2]l. 

During the two orbits the accretion rate onto the two black 
holes is highly variable. The secondary black hole starts ac- 
creting soon after the plunge and is followed by the primary 
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FIG. 6. — Cooling rate of the gas as a function of time calculated from the 
BB-model. The cooling rate shown here traces the bolometric light curve of 
the gas disk (the power generated locally in the disk). The arrows mark the 
times of pericentric passages of the binary. 

which starts accreting a year later at a much higher, super- 
Eddington rate M ~ 10 3 M yr~' (Pig EJ. The jump in accre- 
tion rate is caused by a sudden increase in outward angular 
momentum transport in the disk triggered by a gravitational 
perturbation of the disk by the secondary. We observe such a 
high accretion rate in this simulation only. This happens be- 
cause the high density gas in the central portion of the disk in 
the BB cooling model initially equilibrates at lower tempera- 
tures than in the photoionization models. The cold gas with a 
very low velocity dispersion is more easily driven towards the 
primary black hole when the secondary reaches the pericenter, 
at the end of year 7 of the simulation. The second pericentric 
passage of the secondary occurs 23 years after the beginning 
of the simulation. During the second passage there is no large 
increase in the accretion rate because the gas has acquired a 
significant amount of internal energy and additional kinetic 
energy and it is not easily swept up by the secondary. The 
accretion rates onto the primary and secondary black holes 
become comparable later in the simulation. Unlike the sec- 
ondary, the primary black hole exhibits large variations in the 
accretion rate (^2 orders of magnitude). This occurs because 
the primary black hole accretes colder gas, confined to the 
plane of the disk, while the secondary black hole accretes 
from the hotter, more uniform halo as it sweeps through it, in 
a process resembling Bondi-Hoyle accretion. If the accretion 
continues at the mean rate derived from the first two orbits the 
disk mass would be exhausted in only ^300 years. In order 
for accretion to persist the accretion disk needs to be continu- 
ally replenished by an external supply of gas. 

The cooling rate of the gas closely follows the heating rate 
(Fig [6j, corresponding to a bolometric luminosity of the gas 
of ~ 10 43 ergs -1 . The bolometric luminosity as a function 
of time does not exhibit large fluctuations and remains at a 
nearly constant value after the first 10 years of binary evolu- 
tion. In the BB-model the luminosity of the accretion powered 
sources dominates over the bolometric luminosity of the gas 
disk (i.e., the power locally generated in the disk) by at least 
one order of magnitude during the period of "uniform" accre- 
tion and by several orders of magnitude during the first peri- 
centric passage. This is a consequence of inefficient radiative 
cooling (compared to photoionization calculations). In the BB 
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FIG. 7. — Sequence of snapshots from the simulation showing the evolution of the binary and gas in the S-model (projected into the plane of the binary orbit). 
The time stamps are marked on the top of each panel and are the same as in Figure|4] The first snapshot shows the system shortly after the secondary plunges into 
the disk for the first time. Note that the second pericentric passage in the S-model occurs about 100 days later than in the BB-model. The rotation of the binary 
and the disk is counter-clockwise. The color legend is the same as in previous figures. The higher temperature particles are plotted over the lower temperature 
ones, with the result that some information is hidden (see Figure[8j. 



model, the inefficient radiative cooling is compensated by adi- 
abatic expansion of the gas, resulting in a more pronounced 
spatial dispersion of the gas. 

3.2. Co-rotating Case with Solar Metallicity Gas (S-Model) 

Figure [7] shows the evolution of the morphology and tem- 
perature of gas at the same times as Figure |4] The initial 
steady state disk in the solar metallicity model (hereafter, the 
S-model) has a temperature of ~ 10 3 K and does not exhibit 
the "two phase" temperature distribution, seen in the BB- 
model. The speed of sound in the unperturbed disk is about 
5 kms -1 . The bolometric luminosity of the gas disk during 
this period is powered by collisional excitation only and it is 
of order 10 40 - 10 41 ergs- 1 (FigureH). 

The gas in the S-model shows the same main morpholog- 
ical features as in the BB-model. The differences between 
the two are that after the interaction with the secondary the 
gas in the S-model remains cooler and shows less dispersion 
with the result that the spiral arm and filaments are more pro- 



nounced. Because of the way we plot particles in Figure [7] 
(higher temperature particles are plotted over the lower tem- 
perature ones), part of the information remains hidden. For 
a better illustration, in Figure [8] we plot the temperature dis- 
tribution in the disk as a function of radius which shows that 
multiple temperature components can be present at a single 
radius. The disk density at the end of the simulation is in 
the range 10 10 — 10 14 cm~ 3 , and the median density is about 
one order of magnitude higher than in the BB-model (see Ta- 
ble[2jl. The disk gas in this model also remains well above the 
Toomre threshold for gravitational instability. The speed of 
sound is about 5kms _1 in the cold disk and ~ 10 2 kms _1 in 
the hotter gas component. The temperature of the gas reaches 
the highest value (T~ 10 12 K) after the shock is formed by the 
secondary. On a time scale of months, the temperature of the 
gas falls below 10 10 K due to the combined effects of radia- 
tive cooling and adiabatic expansion. The dominant radiative 
cooling mechanisms of the shocked gas are inverse Comp- 
ton and free-free emission. The hot component of the gas 



12 



Bogdanovic et al. 



Fraction of Particles 
0.0 0.1 0.2 0.3 0.4 



10 1 



io 10 - 



— I 1 1— 





2000 4000 
Radius (r^) 



6000 



FIG. 8. — Left: Temperature of gas cells as a funcion in the disk which was partially hidden in the xv'-projection in Figure [7] The colors that mark each 
temperature band are used throughout the paper. This figure corresponds to the morphology of the disk plotted in panel 2 of Figure^] Right: Histogram showing 
the fraction of particles as a function of temperature. 



spends a significant amount of time in the temperature range 
10 4 - 10 8 K. In this regime the radiative cooling is dominated 
by free-free emission and recombination radiation, while in- 
verse Compton emission becomes inefficient in comparison. 
On the other hand, the cold component of the gas, confined to 
the spiral arm, has the same physical properties as that in the 
BB -model and it retains a temperature of 10 3 — 10 4 K. 

In the period before the impact of the secondary, the radia- 
tive cooling rate of the solar metallicity gas matches that of 
the BB-model (Figs|6]and|9]l. This is an important test of the 
S-model, because steady, geometrically thin disks are indeed 
expected to be optically thick and behave as a black-body. The 
differences in the cooling processes between the S-model and 
the BB-model appear after the (accretion-powered) sources 
of ionizing radiation turn on. The gas in the S-model exhibits 
less efficient heating and more efficient cooling relative to the 



BB-model. This happens because the ionized gas remains op- 
tically thin to ionizing radiation and absorbs it less efficiently, 
while at the same time free-free and recombination cooling 
processes are boosted. Once the gas has been ionized its phys- 
ical properties depart from these in the black-body model. 

After the first pericentric passage of the secondary 
the accretion rate reaches 10M Q yr _1 (Figure ITOb. com- 
parable to the Eddington rate for this system, Me ~ 
30M 8 (tj/0.01) -1 M yr -1 . During the remainder of the sim- 
ulation the accretion rate remains at a nearly constant level 
of about lM yr _1 . This implies that the accretion lumi- 
nosity is just below the Eddington luminosity (Le = 1.51 x 
10 46 Mgergs _1 ) for a few years after the first pericentric pas- 
sage of the binary and later settles at luminosity ~ 10 Lg. 
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FIG. 9. — Cooling rate of the gas as a function of time calculated from the 
S-model. The cooling rate shown here traces the bolometric light curve of the 
gas disk. The arrows mark the times of pericentric passages of the binary. 



FIG. 10. — Effective accretion rate on the primary (solid, red line) and 
secondary (dashed, blue line) black holes calculated from the S-model. The 
accretion rate curves can be translated into UV/X-ray light curves by assum- 
ing that I M . u ' ~ 10 43 ergs~' of UV/X-ray luminosity. The arrows mark 
the times of pericentric passages of the binary. 
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FIG. 1 1 . — Sequence of snapshots from the simulation showing the evolution of the binary and gas in the SR-model (projected into the plane of the binary 
orbit). The time stamp is marked on the top of each panel. The first snapshot shows the system shortly after the secondary plunges into the disk for the first time. 
The rotation of the binary is counter-clockwise and that of the disk i s cl ockwise. Note that higher temperature particles are plotted over the lower temperature 
ones, with the result that some information is hidden (also see Figure [TSV 



During this period the bolometric luminosity of the disk, pow- 
ered by shocks and illumination, is found to be ~ 10 45 ergs _1 
on average, with fluctuations of up to 2 orders of magnitude. 
The bolometric luminosity of the gas disk is comparable to 
that of the photoionization sources. Additionally, the pericen- 
tric passages, in the 7th and 23rd years of the simulation, can 
be easily discerned in the cooling curve (Figure [9). These are 
not noticeable in the cooling curve of the BB-model (Figure[6]i 
because its bolometric luminosity is dominated by emission 
from the diffuse gas component. A common property of the 
BB and S-models is large fluctuations in the accretion rate 
of the primary, a consequence of the accretion of cold gas. 
Moreover, the largest dips in the accretion rate of the primary 
coincide with the apocentric passage, when the gas is tidally 
pulled away from the primary and the primary exhibits the 
slowest orbital motion. 



3.3. Counter- rotating Case With Solar Metallicity Gas 
(SR-Model) 

In this scenario the binary and gas disk rotate in opposite di- 
rections with respect to each other. In FigureQT]we show the 
morphology and the temperature distribution, while in Fig- 
ure [12] we plot the temperature stratification in the disk for 
this model at a single time. During the pre-plunge phase the 
physical properties of the disk are comparable with those in 
the co-rotating S -model. The plunge of the secondary into the 
disk in the counter-rotating solar metallicity model (hereafter, 
the SR-model) is more energetic than in the co-rotating mod- 
els, because the relative velocity of the encounter is higher. 
However, we observe a lower median gas temperature in the 
SR-model relative to the S-model. This apparent discrepancy 
can be explained by the fact that the counter-rotating gas disk 
efficiently shears the shock formed by the impact of the sec- 
ondary black hole. The hot gas from the shock front is quickly 
transported away from the secondary and its thermal energy 
is radiated. When compared to each other, the three models 
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FIG. 12. — Left: Temperature of gas cells as a function of ra dius at 18.1 years in SR-model. This figure shows the extent of temperature stratification in the disk 
which was partially hidden in the ry-projection in FigurefTj] The colors that mark each temperature band are used throughout the paper. This figure corresponds 
to the morphology of the disk plotted in panel 3 of Figure TTll Right. Histogram showing the fraction of particles as a function of temperature. 



exhibit an interesting property: while the median gas temper- 
ature is 3 x 10 5 , 5 x 10 5 , and 10 6 K, the mean temperature is 
10 11 , 10 10 , and 10 9 K, in the SR, S, and BB-models respec- 
tively. 

Another property of the SR-model is that its cooling curve 
does not exhibit peaks at the times of the pericentric passages 
of the binary (Figure[T3l). The peaks noticeable in the cooling 
curve occur, respectively, one and two years after the peri- 
centric passages and are not associated with the peaks in the 
accretion rate (Figure [Pfl). This implies that along with ac- 
cretion there must be an additional mechanism that powers 
the bolometric emission. This second heating mechanism, 
which competes with photoionization, arises from shocks in 
the gas. The bolometric luminosity of the gas disk settles 
at ~ 10 45 ergs -1 and the luminosity of the photoionization 
sources is ~ 10 44 ergs -1 . 

In the counter-rotating model, the gas is driven away from 
the center of the disk, as it acquires more orbital angular mo- 
mentum from the binary, and the hollow region around the 
primary widens with time. This is why during a pericen- 
tric passage there is no characteristic jump in the accretion 
rate of the primary, and accretion occurs by smooth diffusion 
of particles from the disk. However, the accretion rate onto 
the secondary shows peaks corresponding to pericentric pas- 
sages, and smooth undulation in between. Although over the 
course of the two orbits we do not notice a decrease in the 
accretion rate, if the outbound migration of gas continues and 
the low density gap is formed at the center, the accretion on 
the primary is likely to dwindle. Observationally, this implies 
that in these systems, only one of the sources (the secondary) 
may remain active, until its orbit gradually becomes smaller 
than the hollow region in the disk, at which point both black 
holes may stop accreting. Whether both, one, or neither of the 
two sources are seen in emission depends on the evolution- 
ary stage of the binary and the thermodynamic properties and 
supply of the low angular momentum gas in the host galaxy. 

A noticeable feature in the accretion curve of the primary is 
a dip between the 8th and 12th years. This corresponds to the 



time when a large portion of the gas is unbound from the disk, 
as shown in the second panel of FigureQT] The detached part 
of the disk continues to revolve around the primary until it is 
re-captured and wound around. After it has completed several 
orbits, this gas creates a complex set of rings and filaments. A 
fraction of this gas will fall into the primary black hole and 
boost its accretion rate back to ~ 1 M yr -1 . 

3.4. X-ray Light Curves 

We calculate the X-ray luminosity powered by the two dif- 
ferent mechanisms: accretion and bremsstrahlung emission 
from the hot gas. The bremsstrahlung emitting gas is heated 
by both, shocks and photoionization. The accretion rates onto 
the primary and secondary black holes estimated from the 
S and SR-models result in a UV/X-ray luminosity of about 
10 43 ergs -1 and up to 10 44 ergs -1 during times of high ac- 
cretion rate (the detai ls of the a ccretion curves in each model 
are discussed in § !3.2l and § !3.3l and shown Figures[T0land[l4l). 

In this section we describe how the estimate of the con- 
current bremsstrahlung X-ray luminosity was obtained. We 
calculate the thermal bremsstrahlung power per unit volume 
from the population of relativistic electrons hotter than 10 K 
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where Z is the atomic number, n e and n, are the number den- 
sities of electrons and ions, respectively, ~gb is a frequency av- 
erage of the velocity averaged Gaunt factor and is of order 1 . 
The electrons are heated in the process of photoionization of 
the gas by the accretion powered sources and in shocks tak- 
ing place in the gas perturbed by the secondary black hole. 
In order to calculate the resulting bremsstrahlung X-ray lu- 
minosity, we also need to know the volume occupied by the 
gas giving rise to this emission. The calculation of the size 
of the emitting volume from a SPH simulation is compli- 
cated by the finite spatial resolution, as set by the charac- 
teristic size assigned to each gas particle (i.e., the smooth- 
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FIG. 13. — Cooling rate of the gas as a function of time calculated from the 
SR-model. The cooling rate shown here traces the bolometric light curve of 
the gas disk. The arrows mark the times of pericentric passages of the binary. 

ing length, h sm i). In our calculations h sm \ > 10 r g rs 10 14 cm. 
As a result of "smoothing", the bremsstrahlung X-ray emis- 
sion from shocks, which typically occur on scales smaller 
than h m i in our simulations, is not resolved. Consequently, 
the size of the emission region in this case is overestimated, 
while the density and the temperature in the shock are under- 
estimated. In order to bypass this uncertainty we calculate 
the values of some physical properties of the X-ray emitting 
gas following the ana l ytical treatment of shock structure in 
iHollenbach & McKed dl979h . 

First, we attempt to discriminate between the gas particles 
in the photoionized and shocked regions. We assume that 
gas particles that have a Mach number ^# > 10 can par- 
ticipate in shocks and thus contribute to the shock-powered 
bremsstrahlung X-ray emission. From the jump conditions 
we calculate the compression and temperature ratios in the 
shocked gas cell and from these we estimate the post-shock 
density and temperature of the gas. Using the post-shock 
values of the physical parameters we calculate the cooling 
rate (equation ffT8l ) and cooling time for the shocked gas 
cell under the assumption that its dominant cooling mecha- 
nism is bremsstrahlung. The corresponding cooling column 
of gas is then N coo i = riQV s t coo i, where «o is the pre-shock 
density of the gas, v s is the velocity of the shock, and t coo i 
is the cooling time. The values of the cooling column and 
the post-shock density of the hot gas allow us to estimate the 
characteristic length scale of the cooling column, l am i- F° r 
typical values of the parameters in the shocked gas cells in 
our simulations, where n e w n,- ~ n ~ 10 12 — 10 cm~ 3 and 
T ~ 10 8 K, we obtain l coo i ~ 10 cm. We use this informa- 
tion to constrain the X-ray emitting volume within a shocked 
gas cell. Given that 60 — 70% of gas particles fulfill the con- 
dition Jt > 10 and post-shock T > 10 7 K, we find that the 
resulting bremsstrahlung X-ray luminosity from the shocked 
gas is in the range 10 — 10 ergs -1 , with peaks reaching 
10 43 ergs -1 . 

In the second step we estimate the bremsstrahlung X-ray 
luminosity emitted from the photoionized portion of the gas. 
We identify the gas cells responsible for most of the X-ray 
emission powered by photoionization as low density cells with 
temperatures higher than T > 10 7 K. The hot, optically thin, 
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FIG. 14. — Effective accretion rate on the primary (solid, red line) and 
secondary (dashed, blue line) black holes calculated from the SR-model. The 
accretion rate curves in our calculations can be translated into UV/X-ray light 
curves from the emission sources, where I \l . \ r ~ lO^ergs" 1 of UV/X- 
ray luminosity. The arrows mark the times of pericentric passages of the 
binary. 

photoionized component of the gas forms a halo with a me- 
dian density of n ~ 10 7 — 10 8 cm~ 3 . The photoionized re- 
gions are more spatially extended than the shocked regions 
and their characteristic scale is £ a >oi ~ 

I0 3 r g w 10 16 cm. The 
number of photoionized gas particles that contribute to the 
X-ray bremsstrahlung emission at any given time is about a 
few percent of all gas particles. However, the large volume of 
this halo makes up for its small mass and low number density. 
Thus, the resulting bremsstrahlung luminosity from the halo 
is 10 41 — 10 42 ergs~ ! , comparable to the estimated contribu- 
tion of the shocked regions. 

We emphasize that the estimates of the bremsstrahlung X- 
ray luminosity presented here are based on astrophysically 
motivated but simplified assumptions. For example, while 
./M > 10 may be a necessary condition for a gas particle to 
participate in shocks, not every such particle will do so. The 
shocks may be avoided in ordered, laminar flows where the 
trajectories of particles do not intersect. However, we ex- 
pect that the gravitational perturbation from the secondary 
black hole will cause shocks, especially in the inner portion of 
the disk where the dynamical time scale of particles is much 
shorter than the period of the binary and the density of the 
gas is higher. In the case of photoionization, not all hot, low 
density gas particles will be entirely photoionized, making the 
effective emitting volume (or equivalently the emitting mass) 
only a fraction of the total in a gas cell. There may also be 
particles which participate in shocks and are at the same time 
photoionized by the accretion powered sources. For all these 
reasons and because the structure of the nuclear region is very 
complex, we note that the above estimates should be regarded 
as constraints rather than exact values. 

The total bremsstrahlung X-ray luminosity (shock plus 
photoionization powered) of the gas in our simulations 
is lower than the accretion-powered UV/X-ray luminos- 
ity, although during the times of pericentric passage the 
peak bremsstrahlung X-ray luminosity is comparable to the 
accretion-powered X-ray luminosity. The latter X-ray light 
curve exhibits peaks during pericentric passages of the bi- 
nary, thus the peaks are good markers of such events. The 
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FIG. 15. — Ha light curve for the S-model obtained with help of photoion- 
ization calculations. The arrows mark the times of p erice ntric passages of the 
binary. The details of the method are described in i; l2.3l of the text. 

estimated level of the total X-ray emission (accretion pow- 
ered plus bremsstrahlung) should be observable to a redshift 
of z < 2 during the outburst phases. 

Although our calculations follow the evolution of the bi- 
nary over only two orbital passages we suggest that in binary 
systems of this type, X-ray outbursts should be expected to 
continue as long as there is gas in the nuclear region for the 
binary to interact with. Our models imply that at the mean 
level of accretion, the gaseous disks should be depleted after 
10 3 — 10 4 yr. In order for the outburst activity to last for longer 
time scales, the reservoir of gas in the nuclear region needs to 
be continually replenished. It is also plausible, however, that 
repeated collisions of the secondary black hole with the disk 
will completely disrupt it and turn it into a spherical halo of 
hot gas. In such a case, the accretion rate would be relatively 
smooth and uniform, and no outburst would be evident. 

A calculation following the X-ray light curve variability 
over a large number of orbits is necessary in order to confirm 
that the periodicity is a long lived signature of the binary. Al- 
though currently not possible with the SPH method used here, 
simulations of the long term evolution of the binary together 
with hydrodynamics and radiative transfer may be achieved 
in the future by a hybrid approach including complementary 
numerical techniques. In particular, they may be realized by 
combining semi-analytical hydrodynamical calculations with 
emphasis on the long term dynamical evolution of gas with 
methods that focus on hydrodynamics and radiative transfer, 
like the one used here. 

3.5. Ha Light Curves and Broad Emission Line Profiles 

By modeling the Ha light curves for the S and SR-models 
we find that after the beginning of accretion the Ha luminos- 
ity gradually reaches 10 — 10 40 erg s _1 in both models (Fig- 
ures Q3] and [16). Sources with such Ha luminosities are ob- 
servable out to the distance of the Virgo Cluster and possibly 
up to a distance of 100 Mpc. During the steady-state phase, 
before the plunge of the secondary into the disk, both models 
are characterized by a low Ha luminosity, Ly a ~ 10 35 ergs -1 . 
This is an artifact of the calculation, a consequence of the 
choice of initial conditions, because in the pre-plunge phase 
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FIG. 16. — Ha light curve obtained with help of photoionization calcula- 
tions used in the SR-model. The first, large spike in the light curve occurs 
after the first pericentric passage of the binary (the arrows mark th e tim es of 
pericentric passages). The details of the method are described in § |2.3| of the 
text. 

the emission from the gas is entirely thermal emission from 
the disk and the effects of photoionization are not taken into 
account. Such neutral and non-illuminated accretion disks are 
unlikely to be encountered in nuclear regions of AGNs, even 
if these are characterized with very low accretion rates. 

Although in our model the Ha light is closely linked with 
accretion, it does not closely follow the variations in the ac- 
cretion rate. Neither the Ha light curve exhibits easily rec- 
ognizable periodic features during the first two orbits of the 
binary. The only occasion when a noticeable jump in the Ha 
luminosity occurs is in the SR-model, after the first pericen- 
tric passage (Figure [16b. The Ha light is predominantly con- 
tributed by a spatially extended component of the gas and is 
less susceptible to shocks by the secondary. The gradual in- 
crease in the Ha luminosity noticeable in both light curves is 
caused by the redistribution of the photoionized gas in space. 
As gas expands to a scale of a parsec, its density and optical 
depth decrease, and the solid angle subtended to the photoion- 
izing source increases. As a result, a larger fraction of ioniz- 
ing photons is reprocessed to Ha photons and the physical 
properties of the gas evolve towards those found in photoion- 
ized nebulae. 

It appears difficult to infer the existence of a binary in the 
nuclear region of a galaxy based on the Ha light curve alone. 
Serendipitous flares in the Ha light curve may occur close 
to pericentric passages but their duration may be fairly short 
(here, <month, the upper limit set by the time resolution in 
the light curve) and consequently hard to observe. The pe- 
riodic outbursts may be more pronounced in cases when the 
binary orbit is inclined with respect to the disk, but such mod- 
els remain to be explored. 

We also calculate the broad Ha emission-line profiles for 
both the S and SR-models in order to better illustrate the kine- 
matics of the gas (see Figures [TTl and [THI). All profiles have 
been calculated under the assumption that the observer is lo- 
cated at a distance d — > °° in the positive xz-plane, at i = 30° 
to the z-axis. Initially, from the unperturbed disk we observe 
double-peaked emission line profiles but they gradually de- 
part from this shape as the perturbation propagates through 
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FIG. 17. — Sequence of the Ha emission-line profiles selected from the S-model. The intrinsic intensity of profiles is plotted against wavelength. The first 4 
profiles in the sequence are multiplied by a factor of 10, so that they can be represented on the same intensity scale with the other profiles. The corresponding 
time from the beginning of the simulation is plotted next to each profile. The inclination of the plane of the disk with respect to the observer is as marked on the 
figure. The vertical dashed line at 6563 A marks the Ha rest frame wavelength. 



the gas. For reasons discussed above, the regular Ha pro- 
files from a quiescent, cold disk, have a low intensity. The 
profiles appear variable on a time scale of months to years, 
both in shape and width. The profiles in both sets show an 
extended, low intensity red wing, most pronounced during 
the periods of increased accretion and after the times of peri- 
centric passages. During such events a number of emitting 
gas particles find themselves deeper in the potential well of 
the binary. The emitted photons that leave the potential well 
are gravitationally redshifted and Doppler boosted, thus con- 
tributing to the extended red wing and the blue shoulder of the 
relativistic Ha profiles. The changes in the profile sequence 
can be better seen in the trailed spectrograms, shown in Fig- 
ures [19] and [20] which are the 2D maps of the Ha intensity 
against observed velocity (or wavelength). The Ha emission 
in the trailed spectrograms from both S and SR-models ap- 
pears redshifted relative to the rest wavelength. In addition, 
the emission in the extended red wing of the profile is eas- 
ily noticeable, especially in the SR-model (Figurel20ii. where 
this effect is most dramatic after the pericentric passages of 
the binary. Since Doppler boosting of the blue side of the line 
and gravitational redshift of the red wing are pronounced dur- 
ing and after the pericentric passages of the binary, they can 
serve as indicators of the orbital period. In the trailed spec- 
trogram of the S-model it is possible to discern the repeating 
behavior in the windows between 7 to 22 years and 22 to 36 
years. The width of the Ha profile increases after the peri- 
centric passages of the system, reflecting the inflow of gas 



towards the primary. Additionally the widening of the profile 
appears asymmetric and shifted towards the red with respect 
to the pre-pericentric sequence of profiles. This shift is a sig- 
nature of the motion of the accretion disk which follows the 
primary on its orbit. The shift is even more pronounced in 
the trailed spectrogram of the SR-model which also exhibits 
narrower emission line profiles than the S-model. This differ- 
ence arises from the difference in kinematics of the innermost 
portion of the accretion disk surrounding the primary. In the 
SR-model the gas at the inner edge of the disk becomes grad- 
ually less bound to the primary black hole over time due to the 
interactions with the binary. We compare the trailed spectro- 
grams with the velocity curves of the two binary components 
projected onto the line of sight (right panels in Figures [191 
andl20b. Because of the nonzero eccentricity of the orbit the 
velocity curves appear skewed. The variations in the profile 
intensities in the S-model correspond to the features in the ve- 
locity curve of the secondary in year 7 and 23, when it moves 
away from the observer with the highest projected velocity. In 
the SR-model the velocity curve of the secondary appears as 
a mirror image of certain features in the trailed spectrogram 
because in this model the gas rotates in the opposite direction 
with respect to the binary. In general, we expect the signature 
of the velocity curves of the two black holes to appear in the 
emission line profiles at pericentric passages, because that is 
when the binary interaction with the gas is strongest. If such 
features could be discerned in the observed Ha emission-line 
profiles, they would signal the presence of a binary. 
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FIG. 18. — Same as previous figure, except for the SR-model. The first 2 profiles in the sequence are multiplied by a factor of 10, so that they can be represented 
on the same intensity scale with the other profiles. 



In practice, one may constrain the period and the mass ratio 
of the binary from the periodicity and projected velocity com- 
ponents of the two black holes, all measured from the Ha 
profile sequence. The determination of the individual black 
hole masses is more challenging, because it requires knowl- 
edge of the inclination of the binary orbit with respect to the 
observer. In general, the inclination of the binary orbit may 
not be the same as that of the accretion disk but in some cases 
it may be possible to determine it relative to the accretion disk. 
If the binary orbit is circular, fewer parameters are needed to 
constrain the properties of the system, and it should be possi- 
ble to determine the individual masses of black holes, given 
a knowledge of the observed velocity components of the two 
black holes at a single point of time, orbital period, and the 
inclination of the system. To constrain the shape and the ori- 
entation of an elliptical binary orbit one would in addition 
need information about the velocity components at percen- 
ter and apocenter, as well as the orbital separation of the bi- 
nary. Before any of these parameters can be measured from 
observations of light curves and the Ha line sequence, it is 
necessary to follow at least a few revolutions of the binary. 

We note several other qualitative features of the Ha pro- 
files, which can shed light on the motion of gas in the nuclear 
region. In both models, the central part of the profile, con- 
tributed by emission from the low velocity gas, becomes more 
pronounced with time. This low velocity core of the profile 
is contributed by filaments of gas which have been expelled 
from the center at the expense of the angular momentum of 
the binary. The same effect re-appears during the second or- 
bital passage when a new low velocity component is added to 



the system. In this process filaments with a similar morphol- 
ogy are formed during each passage, on a slightly different 
length scale, because the older filaments have had time to ex- 
pand, and with a slightly different phase, due to the effect of 
rotation. 

Although it is expected that all AGNs host black holes fed 
by an accretion disk, broad, double or multi-peaked emission- 
line profiles are observed only in a small fraction of objects 
dEracleous & Halpernl 12001 IStrateva eTail 120031) . The ap- 
pearance of the Ha line profiles can be affected by physical 
conditions in the nuclear region, other than the phase-space 
distribution of the emitting gas, such as optical depth to Ha 
photons and electron scattering. Figure[T]illustrates that a sig- 
nificant fraction of the Ha light is contributed by the hot ion- 
ized gas located in the inner disk, via recombination. This 
component of the gas is transparent to Ha photons (Figure|2|. 
Nevertheless, the profiles will likely be affected by the high 
Ha opacity in regions of the disk where the gas temperatures 
are just below ~ 10 6 K and the neutral hydrogen column den- 
sity is high (Figure |2TI). 

The second factor, which impacts the shape of the Ha 
emission-line profile, but to a lesser extent the Ha luminosity, 
is electron scattering. As expected, the ionized gas will have 
the highest optical depth to Thomson scattering (Figure l22l). 
We find that the maximum value of Xj in our calculation is 
about 10, observed for the hottest gas cells. About 60% of 
the mass of the gas has optical depth to Thomson scattering 
larger than 1, in our simulations. When Thomson scattering 
is significant, the hot corona with a high thermal velocity dis- 
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FIG. 19. — Trailed spectrogram (left) and black hole velocity curves (right) plotted for the S-model. The trailed spectrogram is a logarithmic gray scale map of 
the Ha intensity against wavelength and observed velocity. Darker shades mark higher intensity. Notice the low relative intensity of profiles before the start of 
accretion. The velocity curve panel shows the orbital velocities of the primary (dashed, red line) and secondary (solid, blue line) black holes projected along the 
line of sight to the observer. The velocity curves are skewed because of the non-zero eccentricity of the orbit. 
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FIG. 20. — Trailed spectrogram (left) and black hole velocity curves (right) plotted for the SR-model. The trailed spectrogram is a logarithmic gray scale map 
of the Ha intensity against wavelength and observed velocity. Darker shades mark higher intensity. Notice the low relative intensity of profiles before the start 
of accretion. The velocity curve panel shows the orbital velocities of the primary (dashed, red line) and secondary (solid, blue line) black holes projected along 
the line of sight to the observer. The velocity curves are skewed because of the non-zero eccentricity of the orbit. 

of the Ha emission line profiles, we are not able to account 
for these effects in our calculations. We used conservative val- 
ues for the model parameters, such as the radiative efficiency 
and covering factor, throughout the paper so that the Ha lumi- 
nosity is not significantly over-estimated. The problem of the 
propagation of photons along the line of sight is numerically 



persion (c s ~ 10 2 — 10 3 km s _1 ) could smooth out some emis- 
sion features in the Ha profiles, although the emission profiles 
should still show up clearly. Such conditions can be met in 
nuclear regions hosting a strong photoionization source. 

Although the effects of optical depth in the gas can play an 
important role in estimates of Ha luminosity and appearance 
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FIG. 21 . — Neutral hydrogen column density of the gas as a function of ra- 
dius, at a time of 9.4 years in the S-model. The colors mark the temperatures 
of gas particles and are the same as in earlier figures. This figure corresponds 
to the morphology of the disk as plotted in panel 2 of Figure [7J 

complex and is likely to be accurately treated with the new 
generation of codes which incorporate ray-tracing techniques 
in the radiative transfer calculations. 

3.6. Evolution of the Binary Orbit 

In this section we describe the changes in the binary orbit 
due to interactions with the gas and compare them with an- 
alytical estimates. The evolution of the binding energy and 
orbital angular momentum of the binary in the BB, S, and 
SR-model is shown in Figures [23] [24] and [25] respectively. 
In all three models the binary system shows a net loss of 
both energy and orbital angular momentum. The mean dis- 
sipation rate of the binding energy, calculated from the sim- 
ulations is dE/dt w4 x 10 46 ergs~' for the BB-model and 
1 x 10 46 ergs~' in the S, and SR-model. As a result the sec- 
ond pericentric passage in the S-model occurs about 100 days 
later than in the BB-model (see Figures |4] and [7] at a time of 
23.2 years). The mean rate of energy dissipation due to emis- 
sion of gravitational radiation is [dE /dt) gw K,l x 10 42 ergs _1 
and thus the evolution of the binary orbit is dominated by in- 
teractions with the gas. Because the rate of gravitational wave 
emission is very similar for all three models we show it only 
once, in the inset of Figure [24] 

In all three models the binary is initially placed on the same 
orbit, hence the differences in evolution of the orbital param- 
eters among the models are entirely due to the kinematic and 
thermodynamic properties of the gas. Recall that in the BB- 
model the gas initially redistributes into two phases with a 
temperature contrast of about an order of magnitude (§ 13. 11 1. 
The cold phase present in the BB-model is missing in the 
S-model, where the secondary black hole interacts with the 
warmer and more tenuous gas and consequently dissipates its 
energy more gradu ally. This is consistent with the findings of 
lEscala et alJ (1200 5) that the binary evolution depends strongly 
on the "dumpiness" of the gas, and that the binary suffers a 
stronger deceleration in gaseous media with weaker pressure 
support. After the first pericentric passage, the cold gas phase 



in the BB-model disappears, due to shock and photoionization 
heating, and E and L evolve in a less dramatic way. 

The rates of dissipation of energy and orbital angular mo- 
mentum in Figures [23] [24] and [25] are functions of the or- 
bital phase of the binary because the binary spends a fraction 
of its time moving through the disk. Therefore, the interac- 
tion of the binary an d the accretion disk occur in two different 
regimes studied by iLin & Papaloizoul dl979allbb . In the first 
regime the disk is confined to the Roche lobe of the primary 
black hole and the secondary is outside the disk and interacts 
wi th it tidally. For a mass rat io q and and a Reynolds number 
M ILin & Papaloizoul d!979al) show that when q > £tfr x , the 
tidal force overtakes the viscous expansion of the disk, and 
the secondary effectively removes the angular momentum at 
the outer edge of the disk by truncating the disk (in our sim- 
ulations q = 0.1 and & ~ 10 5 ). The excess angular momen- 
tum is transported to the binary via tidal torques. In the sec- 
ond regime, the secondary is orbiting withi n the circumbinary 
disk an d truncate the inner edge of the disk (ILin & Papaloizoul 
1 1979b!) through its tidal forces. In this regime the disk acts as 
an effective energy and angular momentum sink and the bi- 
nary should evolve towards a more bound state via resonant 
interactions with the disk. This is qualitatively consistent with 
our findings from the BB and S-models, where the binary is 
co-rotating with the disk and is found to reach the minimum 
binding energy when both black holes are close to the peri- 
center and within the disk. 

An additional torque arises from the fact that the veloc- 
ity vectors of the binary on the elliptical orbit and the gas 
around the primary on circular orbits continuously change 
their respective orientation; thus the torque acting on the bi- 
nary changes accordingly. In the co-rotating case, the sec- 
ondary experiences a torque that pushes it away from the pri- 
mary as it approaches pericenter and a torque that pushes it 
inwards as after it passes pericenter. This effect can be seen 
in Figure [24] the binary evolves from a more bound to a less 
bound energy state as it passes pericenter and the slope of the 
energy curve is negative. In the counter-rotating case, the sec- 
ondary and the disk orbit in the opposite sense and the effect 
of torques is reversed. Figure[25]shows this reversal: the slope 
of the energy curve is positive during the pericentric passage. 

We first compare the rate of dissip ation of or bital energy of 
the binary with that predicted by Ostriker ( 1999) for the case 
where a massive perturber moving supersonically experiences 
the dynamical friction force arising from the wake that formed 
in the ambient medium. This rate can be expressed as 

/dE\ = An{GM p fp /r M \^ 
\dtj df v \r m in J 

(19) 

where M p and v are the mass and velocity of the perturber, 
p is the density of the gaseous medium, and r mux and r„„„ 
are related to the linear size of the medium and the perturber. 
If we assume that r mca is comparable to the radius of the 
disk and r m in to the accretion radius of the black hole, then 
^(r majc / r m i„) ss 5. In our simulations the wake of the hot gas 
trailing the secondary black hole can be seen in panel 1 of 
Figures [4] and[7] 

However, as discussed by Lin & Papalo izoul dl979albl) . 
tidal torques and resonant scattering play an important role in 
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FIG. 22. — Left: Optical depth to Thomson scattering of each gas cell as a function of radius, at a time of 9.4 years in S-model. This figure corresponds to the 
morphology of the disk as plotted in panel 2 of Figure|7] Right: Histogram showing the fraction of particles as a function of Thomson optical depth. 



the interaction of a binary with a shear flow. We also note that 
dynamical friction cannot be the dominant dissipation mech- 
anism in the SR-model. In this model, the disk rotates in the 
opposite sense from the binary and the wake formed by the 
secondary is efficiently she ared away from the second ary. We 
utilize the expression from lLin & Papaloizoul d!979bl) for the 
tidal torque exerted by the binary on the disk in the process of 
gap formation and obtain 
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where D. is the angular velocity of the gas interacting with the 
black hole at radius r, E is the surface density of the gas disk, 
A is the half-size of the gap, and H is the half-thickness of the 
disk. The results of our simulations are consistent with both 
mechanisms being responsible for the dissipation of the bi- 
nary orbital energy. We also note that dynamical friction can- 
not be the dominant dissipation mechanism in the case of the 
binary and gas disk rotating in the opposite sense from each 
other. Assuming that the binary maintains a rate of energy 
dissipation of ~ 10 46 ergs _1 , it would proceed to coalescence 
in only ~ 10 3 yr. Most likely, the rate of inspiral due to in- 
teractions with the gaseous backgroun d will not persist at this 
high r ate over long periods of time. Armitage & Nataraian 
(2002) for example discussed the evolution of a binary on 
a circular orbit, in the gas disk, from separations of 0.1 to 
10~ 3 pc. Extrapolating from their results one obtains that 
from the separation of 0.01 pc, a combination of disk-driven 
migration and inspiral due to gravitational radiation leads to 
coalescence within few x 10 6 years. In order to obtain a robust 
estimate of the orbital decay rate one should follow the evo- 
lution of the binary and gas starting from much larger separa- 
tions. While this task was undertaken by a number of authors 



(Kaza ntzidis etai]l20 05: Escala_etaL||2004, 2005: lDotti etail 
I2006at iMaver et al.| ]2006a) it is still a computational chal- 
lenge to extend the wide dynamic range of the simulations 
to the lowest binary separations, where gravitational radiation 
overtakes the dissipation caused by the gaseous background. 

4. DISCUSSION 

Here we discuss the impact of underlying assumptions in 
our models on the observable electromagnetic signatures and 
the relevance of our results in the context of the larger scale 
simulations. 

Mass of the gas disk. A more massive disk would change the 
gas drag, hence make the binary evolve faster and vice 
versa. This would affect the number of pericentric pas- 
sages of the binary and in turn the number of peaks in 
the cooling rate and X-ray light curves occurring over 
a given period of time. The number of shifts in wave- 
length space noticeable in the Ha trailed spectrogram 
would also change, reflecting the change in the orbital 
period of the binary. This may complicate the predic- 
tion and observational followup of the outbursts, since 
after the first few orbital cycles the uncertainty in the 
orbital period of the binary will be rather large (weeks 
to months; see the example of OJ 287). If a decrease in 
the orbital period of a MBHB candidate can be deter- 
mined from the observed light curves, it would imply 
that the binary is dissipating orbital angular momentum 
via interactions with the gas disk or emission of gravi- 
tational waves. 

The MBHBs likely to be discovered by monitoring of 
light curves and Ha emission line profiles of candi- 
date galaxies would likely have orbital periods of order 
of 10 yr or less, implying binaries with subparsec or- 
bital separations. The dominant dissipation mechanism 
could be determined if a constraint on the masses of 
black holes can be put in addition to the rate of change 
of the orbital period. One could then compare the ob- 
served rate of inspiral to that predicted by the theory of 
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FIG. 23. — Top: Evolution of the binding energy of the binary over two 
orbits in the BB-model. Inset: Binding energy in the period 10 to 40 yr 
in more detail. Bottom: Evolution of the orbital angular momentum of the 
binary. The arrows mark the times of pericentric passages of the binary. Inset: 
Evolution of the orbital separation of the binary. 

general relativity for some assumed eccentricity; if the 
two are discrepant this could be explained by presence 
of the gaseous dynamical friction. 

Structure of the disk. Besides the total mass, the efficiency of 
dynamical friction depends on the structure of the ac- 
cretion disk. The intensity of the gas drag is stronger 
in systems where the density of the gas is higher and 
where a massive body moves su personically w ith re- 
spect to the ambient medium (Ostriker 1999). The 
structure of nuclear accretion disks on a subparsec scale 
is rather uncertain and in our calculations we assumed 
that the disk is geometrically thin. The r e sults from 
the larger scale simu lations by lEscala et al.l d2005l) and 
May er etakl ((2006a) provide useful insights into the ef- 
fect that a different equation of state of the gas disk can 
have on the orbital decay of the binary at separations 
of 10—100 pc. Both groups find that the binary inspi- 
ral is more efficient in galactic disks with less pressure 
support and smaller scale heights. Such disks resem- 
ble clump y and dense disks in starforming galaxies. In 
addition, Maveretal. (2006a) find that for a galactic 
merger to result in a bound pair of massive black holes 
and eventually a coalescence within the Hubble time, 
the AGN should not play an important role as an source 
of heat in the disk during the merger phase. 



Role of starburst and AGN. A number of authors found 
that when a large fraction of the disk gas is 
funneled to the central 1 kpc, following a galac- 
tic merger, it can lead to n uclear starburst on 
times cale less than 10 8 yr ( Barnes & Hernquisil 
19961: iDi Matteo. Springe!. & Hernquist 120051: 



Springel, Di Matteo, & Hernquist 2005). Because 



it may take 10 7 yr or longer for a binary to sink from 
kiloparsec to subparsec scales, the MBHB inspiral and 
coalescence in such galaxies may coincide with an 
active starburst dDotti et a l. 2006b). It appears that if 
a starburst reaches its peak before the MBHB merges, 
a large fraction of the gas may be already turned into 
stars while the remaining gas will be exposed to the 
starburst feedback, thus diminishing the role of gas 
as a dissipation agent in the MBHB inspiral. The 
resulting accretion rate onto the black holes may be too 
low to give rise to observable signatures in the form 
of periodic outbursts and variable Ha emission line 
profiles; these galaxies would app ear as ultraluminous 
or luminous s tarburst galaxies (Fa bian et al.l 119981 : 
iDotti e t al. 2006b). However, if the onset of a strong 
AGN activity in either of the parent nuclei precedes the 
starburst, the mechanical and radiative feedback could 
inhibit both the starburst and further accretion onto 
the black holes and thus suppress the electromagnetic 
signatures of the binary. 

Model of radiation physics. Since a fully 3D radiative trans- 
fer calculation has not been developed yet the model of 
radiation physics that is employed here relies on sev- 
eral assumptions. Here we address the key assumptions 
that can directly affect the balance between heating and 
cooling in the simulation, and consequently the elec- 
tromagnetic signatures of the binary system. The first 
assumption is related to the covering factor, £, used to 
parametrize the effects of self-shielding and geometry 
when gas is illuminated by the incident radiation (equa- 
tion [13}. In our calculation this coefficient has a fixed 
value while realistically it is a function of time and the 
position of a particle with respect to the source of il- 
lumination. Varying £ would result in a different tem- 
perature in the gas hence, different amounts of hot and 
cold components of gas. For example, a larger value 
of £ would amount in an increase of the diffuse, hot 
component of gas. This would boost the Ha and X-ray 
luminosities but may not affect significantly the conclu- 
sions about the variability of these components. Alter- 
natively, a smaller value of £ would result in more gas 
being in the cold phase and confined to the plane of the 
disk. In this scenario both the Ha and X-ray luminosi- 
ties would be lower, and the Ha light could also exhibit 
a periodicity. This is because the colder and denser gas 
is more susceptible to shocks induced by the secondary 
black hole. 

Another important assumption in our calculation is that 
radiative cooling of the gaseous accretion disk can be 
approximated by a sum over the discrete gas particles. 
Each particle represents a gas cell of mass 0.1 Mq. 
The approximation should converge to an exact solu- 
tion when Npart —> °° m case when the gas is optically 
thin to radiation. We thus conclude that for the dif- 
fuse gas component in our simulations this approxima- 
tion should be a satisfactory description and we turn 
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FIG. 24. — Similar to previous figure, except for the S-model. Top in- 
set: Energy loss rate due to emission of gravitational waves over time. The 
binary loses orbital energy via gravitational waves at a mean rate of about 
7 x 10 42 ergs~' and emits with the highest rate close to the pericenter. Bot- 
tom inset: Evolution of the orbital separatio n of the binary. 

attention to the optically thick component. This com- 
ponent has a higher neutral hydrogen column density 
(Figure |2"TTi and is confined to a disk with the thickness 
corresponding to one smoothing length, h sm [. There- 
fore, the vertical structure of the disk is not resolved 
and every gas cell within the disk should cool by emis- 
sion from the top and bottom of its surface. This is 
equivalent to a geometrically thin and optically thick 
accretion disk model where gas cools by emission from 
the two faces of the disk. In our model the gas cells 
radiate isotropically, over their entire surface area and 
according to the cooling rate per unit area determined 
by Cloudy. Thus we overestimate the total cooling from 
the optically thick component of gas by a factor of ap- 
proximately 2. The largest error is expected for the BB- 
model, where all gas particles are assumed to be opti- 
cally thick, while the S and SR models should provide 
a more accurate description. 

Role of resolution. Increasing the resolution yields a better 
accuracy in the calculation of angular momentum trans- 
port and at the same time it may cause larger errors 
in the treatment of cooling in our models. We chose 
100k particles for the BB, S, and SR models and car- 
ried out several tests in order to assess the role of reso- 
lution in the heating, cooling, and evolution of angular 
momentum in our simulations. The first set of tests is 
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FIG. 25. — Similar to previous figure, except for the SR-model. 

carried out by evolving the initial conditions (the IC- 
runs, see Table [2j, with the primary black hole and the 
solar metallicity gas disk, before the secondary black 
hole is introduced. Thus we can isolate the dependence 
of angular momentum transport and cooling on resolu- 
tion, in absence of shocks. The results of the three runs, 
IC1, IC2, and IC3 with 20k, 100k, and 500k particles, 
respectively, are shown in Figures [26] and [27] after the 
initial conditions were evolved for 8 years. The evolu- 
tion of angular momentum appears strongest in the 20k 
run where the pressure gradient drives the evolution of 
the surface density at the inner disk radius (see the top 
panel of Figure |26*1 > while at larger radii the distribu- 
tions appear consistent. The dominant component of 
pressure is thermal gas pressure; the viscous pressure 
is of lesser importance and it operates on much longer 
time scales. The accretion disks at lower resolution thus 
appear to have more gas pressure support than at higher 
resolution. Lower resolution also results in more vis- 
cous disks and shorter viscous time scales, as shown in 
middle panel of Figure [26] where we plot the viscous 
time scale, T viscous ~ r 2 /(a Gadget h sml c s ). 

The bottom panel of Figure [26] shows the temperature 
profile of the disk calculated at three different resolu- 
tions, from the same snapshot in time. The disk in 
the 500k IC3 is the coldest compared to the 100k IC2 
and 20k IC1 runs. This is a consequence of the de- 
pendence of the cooling model on the resolution. In 
the solar metallicity models, the disk cooling rate is 
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FIG. 26. — The effect of resolution on the disk surface density (top), the 
viscous time scale (middle), and the temperature profiles (bottom) in models, 
where the secondary black hole was not introduced to the simulation. The 
profiles shown are for runs 20k IC1 (dash-dot,blue line), 100k IC2 (solid, 
black line), and 500k IC3 (dashed, red line), all at a time of 8 years after the 
beginning of a simulation. The gas was assumed to have a solar metallicity. 
See discussion in §|4]for more details. 

expected to scale with the total number of gas parti- 
cles and their emitting volume as °= N part h^ ml — const. 
Figure [27] shows the cooling curves of an unperturbed, 
passively cooling disk for the three runs with the cool- 
ing being strongest in the 500k one. The scaling no- 
tably departs from the analytical expectation because 
the adaptive smoothing length, h sm i, does not decrease 



with the resolution as quickly as ~ N, 
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FIG. 27. — Dependence of the cooling rate on resolution in a model where 
the secondary black hole is absent and the solar metallicity gas disk is cooling 
passively. The curves are from runs 20k IC1 (dash-dot, blue line), 100k IC2 
(solid, black line), and 500k IC3 (dashed, red line), at the time 8 years after 
the beginning of each simulation. 
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FIG. 28. — Effective accretion rate on the primary (solid, red line) and 
secondary (dashed, blue line) black holes calculated from th e 20 k S 1 run, a 
lower resolution equivalent of the 100k S-model run (Figure [Tol . Note that 
this model has been run for a longer time (close to 3 orbital periods) and that 
three peaks are noticeable in the accretion curve. The arrows mark the times 
of pericentric passages of the binary. 



lution ones. Because the cooling rate of a black-body 
emitter is characterized by its emitting surface area, 
rather than volume, the dependence of the cooling on 
resolution is different in the BB-model. In this case 
the cooling rate is nominally expected to increase with 



resolution as °c N par t h] ml 



quently, the higher resolution runs are characterized by 
a larger cooling volume with respect to the lower reso- 



1 /3 

•'Npart ( see S 12.3b . However, 
the cooling rate may have stronger dependence on N par t 
due to the behavior of h sm i with resolution, as described 
previously for the solar metallicity case. 

In order to estimate how much the cooling rates mea- 
sured from our simulations suffer from the systemat- 
ics discussed in this section, we compare the cooling 
rates for the passively cooling disk from the IC runs 
with the analytically calculated cooling rate for a black 
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FIG. 29. — Sequence of Ha emission line profiles c alcu lated from the 20k SI run, a lower resolution equivalent of the 100k S-model run. The emission line 
profiles are plotted for the same time slices as in FigurefTT] for comparison. Note the lower intensity scale for the profiles in 20k S 1 run but otherwise very similar 
profile widths and shapes to those in 100k run. The inclination of the plane of the disk with respect to the observer is as marked on the figure. The vertical dashed 
line at 6563 A marks the Ha rest frame wavelength. 

body with T = 1000K and surface area equal to that 
of the disk. The comparison is meaningful because the 
cold, nearly isothermal, optically thick and geometri- 
cally thin disk can be described as a black body until 
it is perturbed or illuminated. We obtain a cooling rate 
of Lbb = 3 x 10 41 ergs -1 , which is slightly higher but 
within an order of magnitude of the cooling rate in the 
100k IC2 run (Figurel27li. The same is true for the early 
portions of the cooling curves in the BB, S, and SR 
models, before the passage of the secondary through 
the disk (Figures [6] [9] and[T3l. Their difference from 
the analytic estimate value can be explained by slight 
departures of the temperature of the gas particles from 
1000K. We conclude that the cooling rates calculated 
in our models with 100k particles are in fair agreement 
with the analytic estimate in the limit of an optically 
thick gas disk. We also discuss a possibility for fur- 
ther improving the BB cooling model and suppress the 
dependence of cooling on resolution in § 14.11 



In the second set of tests we examine the effect of reso- 
lution on the accretion rates and the variable Ha emis- 
sion line profiles in the presence of the secondary black 
hole. These quantities are connected to the resolution 
through the graininess of the gas and through the ar- 
tificial viscosity, which leads to different results for 
shocks. For the purpose of the test we carried out the S 1 
simulation of the binary with a prograde disk using 20k 
particles (analogous to the 100k S model). The accre- 



tion rate curves (Figure 1281 ) show the same prominent 
features as the ones in the 100k S-model (Figure |T0b- 
The total accretion rate in the 20k S 1 run is in the range 
1-10 M yr -1 and is slightly higher than in the 100k S 
run. In the lower resolution run the contrast between the 
accretion rate of the secondary black hole at pericenter 
and apocenter appears higher, creating pronounced pe- 
riodic humps in the accretion curve. This is likely the 
result of the Bondi-like accretion on the secondary cou- 
pled with a lower particle number density in the 20k S 1 
run. 

In addition, we calculate the Ha emission line profiles 
from the 20k S 1 run in the same way as in the S model; 
the results are shown in Figure [29] As a consequence 
of the lower cooling in the 20k run, the intrinsic inten- 
sity of the Ha lines is several times lower than in the 
100k S-model run but the integrated Ha luminosity is 
the same, L# a ~ 10 40 ergs -1 . This can be understood 
having in mind that while the intrinsic line intensity de- 
creases at lower resolutions, the smoothing radius used 
to estimate the emitting surface of a gas particle, in- 
creases, and the two effects cancel out in calculation of 
luminosity. The shape and width of the profiles show 
the same dominant features as the profiles in 100k S- 
run and do not appear to be very sensitive to the change 
of resolution. 

While employing a higher number of particles has a 
clear advantage we find that a satisfactory accuracy is 
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achieved with 100k particles given the time scales in- 
volved in our simulations. The evolution of angular 
momentum in the 100k runs is slow and the viscous 
evolution even slower and they do not result in shocks 
and transients in the gas. The estimated systematic er- 
ror in cooling rate due to the dependence on resolution 
is within one order of magnitude and the shapes of the 
Ha emission line profiles appear robust. 

Role of accretion model. In the last paragraph of § [2] we de- 
scribe the approximation we used in order to estimate 
the accretion rates of the two black holes given the 
properties of their perturbed but unresolved accretion 
flows. We tested whether the accretion rate onto the 
secondary black hole affects the results by lowering it 
by a factor of 100 in one of the models. We refer to this 
as the S2-model. From this run, in which we employed 
20k particles, we calculated the observational signa- 
tures and found that they are very similar to the 20k S 1 
and 100k S-models (specifically, the cooling curve fea- 
tures and periodicity, Ha profiles, and Ha light curve). 
The implication is that varying the accretion rate effi- 
ciency of the secondary black hole with respect to that 
of the primary has a weak effect on the resulting ob- 
servational signatures. As long as the primary black 
hole maintains equal or higher accretion rate than the 
secondary, the latter can be regarded as a passive per- 
turber. We also note that the estimates of accretion rates 
in both accretion models are rather conservative, thus it 
is more likely that we underestimate rather than overes- 
timate the accretion rate. 

In our simulations we modeled the binary with the semi- 
major axis of 0.015 pc, in the regime when the dissipation due 
to the gaseous background is still dominant over the gravita- 
tional radiation. This regime corresponds to the disk-driven 
migration stage in the 2 D grid simulations of a bina ry and 
locally isothermal disk by lArmitage & Natara ian (2002| here- 
after, AN). Because there is an overlap in the scales treated 
by our respective simulations we compare the results of the 
two works. Note that in our simulations the binary orbit is 
assumed to be elliptical and the disk surface density initially 
corresponds to ~ 10 5 gcm~ 2 at the distance of 100 r g from 
the primary, in comparison to the circular orbit and the sur- 
face density of ~ 10 6 gcm~ 2 in AN. We find that the low den- 
sity quasi-spherical outflow in our models reaches the scale 
of 1 pc in about 40 years, implying the outflow velocities of 
~ 10 4 kms _1 ,in agreement with the prediction of AN. Unlike 
AN, we find that in the phase of the disk-driven migration the 
accretion rate of the primary is higher than the rate expected 
in the absence of a binary. The rate expected in the latter case 
is ~ 10~ 3 M yr _1 while the rate measured from our sim- 
ulations > 10~ 2 M yr _1 . We attribute this difference to a 
temperature distribution of the gas and the eccentric binary 
orbit. In our simulations, in general case, different gas phases 
can be found to reside at the same radii in the disk (Figures [8] 
and[T2"l). The cold component of gas at T^ 10 3 K is mostly 
confined to the plane of the binary orbit where it gets swept 
by the secondary black hole during every passage through the 
disk and accreted. In addition, a binary on an eccentric orbit 
effectively perturbs the accretion disk and drives the spatially 
extended high velocity outflows. The outflows carry away the 
excess angular momentum, allowing the gas at the innermost 
region of the disk to be accreted. 



4.1. Future Prospects 

The physical parameter space of a binary and gas disk sys- 
tem is very large and a thorough parameter study is beyond 
the scope of this paper. In the future it will be useful to inves- 
tigate systems with different mass ratios in which the binary 
orbit and gaseous disk are not coplanar. In these systems the 
secondary black hole does not gradually migrate through the 
plane of the disk but crosses through the disk at some incli- 
nation. Consequently, such encounters are expected to have 
unique emission signatures. The next, even more general set 
of systems to be studied are binaries with two gas disks at 
different relative inclinations. 

Because the model of radiation physics plays a decisive 
role in determining the character of the observational signa- 
tures, we outline a possible avenue for future improvement of 
the BB- model. Introducing a flux-limited diffusion approx- 
imation dClearv & Monaghan|[l999l: IWhitehouse et alj|2005l: 
MayexeLaL||2P06b) can offer a more accurate treatment of 
the optically thick component of gas by emulating a contin- 
uous diffuse radiation field in the inter-particle space. In this 
approximation the SPH particles which are embedded within 
some gas distribution diffuse heat to their neighbor particles 
as determined by the temperature gradient and coefficient of 
diffusion. The particles defining a surface boundary of an 
emitting region in addition to diffusion can freely radiate as 
a black-body or according to some other prescription. The 
main challenge in the implementation is determining which 
gas particles are constituents of the boundary layer in general 
case, when the normal of an emitting surface is not aligned 
with any preferred axis in the system and when it s orientation 
varies in space and time (see Mayer et al. 2006b, for the case 
of a disk geometry). 

Because our simulations are computationally expensive, we 
were not able to follow the longer-term evolution of the binary 
and the disk on time scales on which a gap is expected to form 
in the disk. The gap is a low density region that forms at the 
center of an accretion disk, on a scale comparable to that of 
the binary orbit, as a cons equence of tidal torques that the bi- 
nary exerts on the disk dLin & Papaloizoulll979allbh . While 
we observed a lower peak density in the innermost region of 
the disk at the end of simulations, whether a gap in the disk 
is efficiently cleared or continuously replenished by the gas 
will depend on the gas temperature and hence, on the treat- 
ment of heating and cooling. The simple, exploratory models 
described here can be of importance in cases when the gap 
is either replenished or its formation process is not efficient. 
However in the future we will explore models in which the 
gap has already formed in the disk and investigate its implica- 
tions for the accretion rates and observational signatures. 

5. SUMMARY AND CONCLUSIONS 

We developed a method for calculation of the observational 
signatures of MBHBs with an associated gas component. The 
method comprises SPH calculations carried out with a modi- 
fied version of the code Gadget coupled with calculations of 
cooling and heating of the gas. We implemented and tested 
two different schemes for heating and cooling of gas in which 
the gas is described as a black-body emitter or a photoionized 
gas with solar composition. In the solar metallicity model 
heating and cooling of the gas are pre-calculated with the pho- 
toionization code Cloudy, along with the emission efficiency 
of the Ha line, the optical depth to Thomson scattering, and 
the neutral hydrogen column density. This information allows 
us to calculate the Ha light curves and Ha emission line pro- 
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files from this system. 

Since MBHBs may spend the largest fr action of their lives 
at sub-parsec scales (Begelm an et al.l [l980) it is particularly 
interesting to study the observational signatures of such bi- 
naries. In this paper we investigated the electromagnetic sig- 
natures of binaries based on a few exploratory models; we 
outline the most important findings below: 

• Based on the first set of models, presented here, we find 
that the X-ray outbursts should occur during pericentric 
passages of a coplanar binary. We suggest that period- 
ical X-ray outbursts should persist as long as there is 
a supply of gas that the MBHB can interact with. A 
calculation of a much larger number of binary orbits 
is needed in order to confirm that the periodic outbursts 
are a long lived signature of the binary. At the estimated 
level of X-ray luminosity the binary systems interacting 
with the gas should be observable to a redshift of z < 2. 

• Except for the recurrent outbursts in the X-ray light 
curve the signature of a binary is most easily discernible 
in the Ha emission line profiles. The variable shape 
of the broad Ha emission line profiles can be used 
as a first indicator when searching for MBHBs in the 
nearby universe, within 100 Mpc. Once the binary can- 
didates are selected from a large spectroscopic survey 
(e.g., the SDSS) they can be monitored over long time 
intervals to look for the time-dependent signature of a 
binary. Based on our simulations (for a mass ratio of 
0.1), the wavelength shift in the Ha emission lines over 
the course of an orbital period should be measurable. 
If one can follow the regular variations of the line over 
a few cycles, one could constrain the properties of the 
binary, such as orbital period and the binary mass ratio. 

• The observed signatures predicted here arise from the 
gas disk and accretion onto the two nuclei, both of 
which may not be detectable in galaxies where the nu- 
clear region is obscured or gas exhausted by the star- 
burst. The resulting signatures are in addition sensitive 



to the adopted model of radiation physics but are not 
very sensitive to the number of particles employed in 
the SPH simulations. 

• The evolution of the intermediate phase binaries con- 
sidered here is driven by interactions of the black holes 
with the accretion disk gas, while the emission of grav- 
itational waves is secondary (3-4 orders of magnitude 
lower). Orbital energy is dissipated primarily via dy- 
namical friction against the gaseous background and 
tidal torques in case of a binary co-rotating with the 
disk. In case of a counter-rotating binary, the tidal 
torques play the dominant role in the dissipation since 
dynamical friction becomes ineffective. 

• The current limitation of the numerical method used in 
this study is that the calculation of a large number of or- 
bits following the hydrodynamics and radiative heating 
and cooling becomes prohibitively expensive. There- 
fore, we are only able to model the observational sig- 
nature of the system for a brief period of its evolution, 
after we have assumed a specific initial configuration. 

If the signatures of binaries are found in observations, they 
could be used to estimate the number of MBHBs in this evo- 
lutionary stage and whether MBHBs indeed evolve quickly 
through the last parsec. Although the intermediate phase bi- 
naries themselves emit gravitational waves at frequencies too 
low to be detected with LISA their discovery is important for 
understanding the evolution of the diverse supermassive black 
hole binary population. 
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